MODELING  RELIABILITY  GROWTH  IN  ACCELERATED  STRESS  TESTING 


DISSERTATION 

Jason  K.  Freels 
Major,  USAF 

AFIT-ENS-DS-13-D-02 


DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 


Wright-Patterson  Air  Force  Base,  Ohio 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


AFIT-ENS-DS-13-D-02 


The  views  expressed  in  this  dissertation  are  those  of  the  author  and  do  not  reflect  the 
official  policy  or  position  of  the  United  States  Air  Force,  Department  of  Defense,  or  the 
United  States  Government. 


AFIT-ENS-DS-13-D-02 


MODELING  RELIABILITY  GROWTH  IN  ACCELERATED  STRESS  TESTING 


DISSERTATION 


Presented  to  the  Faculty 

Department  of  Systems  Engineering  and  Management 
Graduate  School  of  Engineering  and  Management 
Air  Force  Institute  of  Technology 
Air  University 

Air  Education  and  Training  Command 
In  Partial  Fulfillment  of  the  Requirements  for  the 
Degree  of  Doctor  of  Philosophy  in  Systems  Engineering 


Jason  K.  Freels 
Major,  USAF 
December,  2013 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


AFIT-ENS-DS-13-D-02 


MODELING  RELIABILITY  GROWTH  IN  ACCELERATED  STRESS  TESTING 


Jason  K.  Freels 
Major,  USAF 


Approved: 

Dr.  Joseph  J.  Pignatiello,  USAF  (Chair]  Date 

Dr.  Raymond  R.  Hill,  USAF  (Member]  Date 

Dr.  Richard  L.  Warr,  Lt  Col,  USAF  (Member]  Date 

Dr.  Joseph  R.  Wirthlin,  Lt  Col,  USAF  (Member]  Date 


Dr.  Adedeji  B.  Badiru  Date 

Dean,  Graduate  School  of  Engineering  and  Management 


AFIT-ENS-DS-13-D-02 


Abstract 

Qualitative  accelerated  test  methods  improve  system  reliability  by  identifying  and 
removing  initial  design  flaws.  However,  schedule  and  cost  constraints  often  preclude 
sufficient  testing  to  generate  a  meaningful  reliability  estimate  from  the  data  obtained  in 
these  tests.  In  this  dissertation  a  modified  accelerated  life  test  is  proposed  to  assess  the 
likelihood  of  attaining  a  reliability  requirement  based  on  tests  of  early  system  prototypes. 
Assuming  each  prototype  contains  an  unknown  number  of  independent  competing  failure 
modes  whose  respective  times  to  occurrence  are  governed  by  a  distinct  Weibull  law,  the 
observed  failure  data  from  this  qualitative  test  are  shown  to  follow  a  poly-Weibull 
distribution. 

However,  using  an  agent-based  Monte  Carlo  simulation,  it  is  shown  that  for  typical 
products  subjected  to  qualitative  testing,  the  failure  observations  result  from  a 
homogenous  subset  of  the  total  number  of  latent  failure  modes  and  the  failure  data  can  be 
adequately  modeled  with  a  Weibull  distribution.  Thus,  the  projected  system  reliability 
after  implementing  corrective  action  to  remove  one  or  more  failure  modes  can  be 
estimated  using  established  quantitative  accelerated  test  data  analysis  methods.  Our 
results  suggest  that  a  significant  cost  and  time  savings  may  be  realized  using  the  proposed 
method  to  signal  the  need  to  reassess  a  product’s  design  or  reallocate  test  resources  to 
avoid  unnecessary  maintenance  or  redesigns.  Further,  the  proposed  approach  allows  a 
significant  reduction  in  the  test  time  and  sample  size  required  to  estimate  the  risk  of 
meeting  a  reliability  requirement  over  current  quantitative  accelerated  life  test  techniques. 
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Additional  contributions  include  a  numerical  and  analytical  procedure  for  obtaining  the 
maximum  likelihood  parameter  estimates  and  observed  Fisher  information  matrix 
components  for  the  generalized  poly-Weibull  distribution.  Using  this  procedure,  we  show 
that  the  poly-Weibull  distribution  outperforms  the  best-fit  modified  Weibull  alternatives  in 
the  literature  with  respect  to  their  fit  of  reference  data  sets  for  which  the  hazard  rate 
functions  are  non-monotone. 
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MODELING  RELIABILITY  GROWTH  IN  ACCELERATED  STRESS  TESTING 


I.  Introduction 

“Our  record  of predicting  where  we  will  use  military  force  since  Vietnam  is 
perfect  -  we  have  never  once  gotten  it  right,  there  isn ’t  a  single  instance: 

Grenada ,  Panama,  the  first  Gulf  War,  the  Balkans,  Haiti,  you  can  just  keep  going 
through  the  list,  where  we  knew  and  planned  for  such  a  conflict  six  months  in 
advance.  ” 

Secretary  of  Defense  Robert  Gates,  May  24,  2011 

Background 

Recapitalizing  an  aging  Air  Force  inventory  requires  a  balanced  approach  in 
which  neither  major  combat  operations  against  near-peer  technology  nor 
operations  such  as  those  employed  in  Afghanistan  and  Iraq  are  over-emphasized 
[1].  With  this  in  mind,  then  Secretary  of  Defense,  Robert  Gates  [2],  unveiled  seven 
strategic  priorities  to  focus  future  science  and  technology  investments  and  guide  the 
development  of  next  generation  of  systems  engineering  tools  and  processes. 
According  to  Neches  [3],  the  Secretary’s  second  priority,  Engineered  Resilient 
Systems  (ERS),  comprises  efforts  to  "efficiently  create,  field  and  evolve  systems 
which  can  readily  adapt  to  the  inevitable  changes  in  threat,  technology  and  mission 
environments.”  Under  the  ERS  construct,  future  systems  must  be  robust  to  a  wide 
range  of  possible  operational  environments,  not  just  a  single  known  scenario.  An 
implication  of  this  strategy  is  the  need  for  critical  components  capable  of  "plug  and 
play"  operations  across  platforms  and  operating  environments. 
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Recent  history,  however,  indicates  that  US  defense  systems  do  not  satisfy  their 
operational  suitability  requirements  sufficiently  often  [4].  Consequently,  the 
Department  of  Defense  spends  too  much  for  system  redesigns,  spares  management, 
and  maintenance.  Reports  commissioned  to  investigate  the  root  causes  of  these 
reliability  and  suitability  shortfalls  [5,  6]  identified  a  lack  of  appropriate  systems 
engineering  processes,  specifically  a  robust  reliability  growth  strategy  beginning 
early  in  the  development  cycle  as  a  primary  contributor.  In  response  to  these 
reports,  the  Reliability  Improvement  Working  Group  adopted  [7,  8]  to  align 
Department  of  Defense  policy  with  the  best  practices  of  reliability  management,  and 
to  provide  the  most  value  with  the  least  risk  in  terms  of  fielding  reliable  products. 
These  standards  specify  a  scientific  approach  to  design  and  build  reliability  into 
products  early-on  and  institutionalize  the  creation  of  a  comprehensive  reliability 
growth  strategy  throughout  the  acquisition  cycle.  As  result  of  their  implementation 
a  greater  burden  is  placed  on  verifying  the  maturity  of  early  designs,  thereby 
minimizing  the  expenditure  of  test  resources  in  subsequent  development  phases. 

Reliability  growth  is  the  positive  improvement  in  a  product’s  reliability 
distribution  parameter  over  a  period  of  time  due  to  changes  in  product  design  or 
manufacturing  processes  [9].  Reliability  growth  modeling  has  historically  played  a 
role  in  determining  whether  major  development  efforts,  such  as  military  weapon 
systems,  are  likely  to  meet  reliability  requirements  in  time  for  graduation  to  the 
next  development  phase,  and  eventually  to  operational  testing.  To  assess  product 
reliability,  prototypes  are  subjected  to  a  series  of  tests  which  exercise  the  system 
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using  a  selected  subset  of  inputs  from  the  overall  set  of  possible  inputs  that  may  be 
encountered  during  the  product’s  lifetime.  Inputs  may  represent  the  expected 
usage  environment  or  comprise  a  specific  stress  profile  directed  by  the  customer. 

Introducing  a  complex  product  often  requires  a  lengthy  development  process 
during  which  it  is  expected  that  the  reliability  will  steadily  improve  based  on  testing, 
failure  mode  discovery,  root-cause  analysis  and  design  changes  or  component 
substitutions.  Testing  may  be  composed  of  many  different  types  of  tests,  each  with 
its  own  objectives.  Developmental  tests  identify  the  technical  capabilities  and 
limitations  of  proposed  designs  and  ensure  sufficient  design  maturity  is  achieved 
prior  to  operational  testing.  In  operational  testing  the  focus  is  on  demonstrating 
that  the  design  is  suitable  for  its  intended  use  in  a  realistic  operational  environment. 
Entering  operational  test  with  an  immature  design  often  results  in  continued 
debugging  into  the  early  life  of  the  product  after  it  has  been  released  to  the  market, 
usually  at  much  greater  cost  than  if  the  fault  were  discovered  in  developmental 
testing.  Therefore  it  is  desirable  to  model  the  improvement  in  reliability  over  time 
to  (1)  forecast  the  length  of  the  development  process,  (2)  ensure  proper  allocation 
of  testing  resources  and  (3)  estimate  the  reliability  upon  entry  into  market. 

Ideally,  the  assessment  of  reliability  growth  should  begin  soon  after  program 
initiation  with  the  development  of  an  idealized  reliability  growth  planning  curve 
(Figure  1).  Once  testing  begins,  progress  can  be  gauged  by  comparing  quantitative 
assessments  of  system  reliability  against  the  idealized  curve.  Close  agreement 
between  the  planned  and  demonstrated  reliability  is  indicative  of  a  successful 
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reliability  growth  program.  Assessments  below  the  idealized  curve  indicate  that 
reliability  is  lagging  and  may  signal  the  need  to  reallocate  test  resources  or 
reevaluate  the  testing  strategy. 


Test  Phase  1  Test  Phase  2  Test  Phase  3 


Figure  1  -  Illustration  of  idealized  reliability  growth  planning  curves  plotted  along 
with  assessed  growth  across  three  test  phases 

Many  defense  acquisition  programs,  however,  fail  to  get  on  their  reliability 
growth  planning  curves  due  to  the  existence  of  too  many  failure  modes  [4]  at  the 
start  of  developmental  testing.  Prior  to  formal  testing  it  is  assumed  that  few  failure 
modes  remain  in  the  system  for  which  the  root  cause  has  not  been  identified  and 
understood.  Mil-HDBK-189  [9]  indicates  that  a  key  predictor  of  program  success 
with  respect  to  meeting  the  reliability  goal  through  growth  testing  is  that  the  initial 
system  reliability  be  at  least  30%  of  the  required  reliability.  Achieving  this  level  of 
initial  reliability  requires  conducting  design  for  reliability  (DfR)  activities  to 
discover  and  remove  both  functional  and  operational  failure  modes  [10]  from  early 
system  prototypes.  Functional  failures  occur  when  a  latent  design  flaw  activates  to 
become  a  failure  during  normal  operation.  Operational  failures  result  when  the 
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operating  environment  to  which  the  system  is  exposed  falls  outside  of  the  expected 
normal  range.  With  the  availability  of  specialized  materials  and  sophisticated 
manufacturing  processes,  functional  failures  have  become  rare  for  all  but  the  most 
complex  systems.  Operational  failures,  however,  are  more  difficult  to  predict  and 
can  result  from  poorly  defined  requirements,  uncertain  operating  conditions  and 
complex  multi-stress  interactions. 

Qualitative  accelerated  stress  testing  (QAST)  was  developed  to  address  both 
functional  and  operational  failures  by  exposing  products  to  elevated  stresses  and 
accelerate  the  identification  and  removal  of  these  failure  modes.  QAST  techniques, 
such  as  [11, 12, 13],  employ  a  test-fix-test  reliability  improvement  strategy  to 
expand  a  product’s  operating  limits  and  thereby  ensure  the  highest  reliability  in  the 
designed  environment  while  potentially  enabling  usage  in  unplanned  environments. 

Problem  Statement,  Objectives  and  Scope 

The  most  well-known  QAST  technique,  highly  accelerated  life  testing  (HALT), 
utilizes  a  step-increasing  stress  profile  to  discover  failure  modes  quickly.  While 
published  case-studies  [14, 15, 16]  demonstrate  the  effectiveness  of  HALT  to 
improve  product  reliability,  no  published  models  exist  to  quantify  the  reliability 
improvement  gained  as  a  result  of  its  use.  Meeker  [17]  questions  whether  sufficient 
data  are  produced  during  HALT  to  construct  a  meaningful  estimate  of  a  system’s 
reliability  improvement.  Typically,  qualitative  tests  are  conducted  on  very  small 
samples  of  early  system  prototypes  (typically  4  -  10)  to  iteratively  find  and  remove 
failure  modes  at  elevated  stresses  resulting  in  successively  more  mature  designs. 


5 


The  product  improvement  process  used  in  HALT  closely  resembles  a  test-fix-test 
reliability  growth  test,  therefore  it  is  hypothesized  that  an  early  estimate  of  system 
reliability  can  be  derived  by  modeling  HALT  as  a  reliability  growth  test. 

In  this  dissertation  a  reliability  growth  projection  model  is  developed  to 
estimate  system  reliability  as  a  result  of  implementing  corrective  actions  to  remove 
failure  modes  exposed  in  accelerated  stress  environments.  The  projection  model  is 
the  first  to  utilize  failure  data  obtained  from  qualitative  accelerated  tests  and 
provides  a  statistically  rigorous  and  defensible  measure  of  the  likelihood  that  a 
complex  system  or  subsystem  can  attain  its  reliability  requirement  within  an 
allocated  test  time.  As  a  result  of  incorporating  failure  data  obtained  at  elevated 
stresses  into  growth  testing,  the  projection  model  enables  a  significant  reduction  in 
the  number  samples  and  the  total  test  time  required  to  estimate  the  reliability  of  the 
improved  system.  This  ensures  that  the  Department  of  Defense  (DoD)  resources 
allocated  for  testing  can  be  used  to  address  problems  that  may  pose  a  significant 
risk  to  system  performance  after  fielding.  As  consequence,  the  cost  of  implementing 
redesigns  for  fielded  systems  may  also  be  avoided.  This  is  especially  important  as 
DoD  resources  are  likely  to  be  reduced  in  the  future. 

Dissertation  Outline  and  Research  Impacts 

This  dissertation  follows  the  scholarly  article  format  where  Chapters  III,  IV  and  V 
represent  stand-alone  papers  that  have  either  already  been  submitted  for 
publication  or  will  be  submitted  after  graduation.  These  Chapters  sequentially 
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develop  the  proposed  model  for  assessing  the  reliability  growth  of  a  product  with 
independent  competing  failure  modes  exposed  to  accelerated  stresses. 

The  main  contribution  of  this  dissertation  is  the  first  model  to  assess  the 
likelihood  of  attaining  a  reliability  requirement  through  reliability  growth  testing 
where  the  test  environment  is  not  representative  of  the  anticipated  field  operating 
conditions.  A  consequence  of  merging  the  accelerated  testing  and  reliability  growth 
methodologies  is  that  one  or  more  of  the  implicit  assumptions  associated  with  each 
individual  methodology  is  violated.  Chapter  III,  therefore,  outlines  several 
accelerated  test  data  analysis  techniques  currently  in  use  and  illustrates  why  these 
techniques  are  insufficient  to  model  data  obtained  from  a  qualitative  accelerated 
stress  test.  A  framework  to  guide  future  research  efforts  is  then  presented  along 
with  specific  next  steps  detailing  ways  in  which  data  from  these  qualitative  tests 
may  be  incorporated  into  a  reliability  estimate  for  a  complex  system. 

In  Chapter  IV,  it  is  shown  that  for  systems  with  independent  competing  failure 
modes  whose  respective  times  to  occurrence  are  each  governed  by  a  distinct 
Weibull  law,  the  observed  system  failure  times  follow  a  poly-Weibull  distribution 
with  vector  valued  parameters  a  and  /?.  A  numerical  and  analytical  procedure  is 
derived  for  obtaining  the  maximum  likelihood  parameter  estimates  and  standard 
errors  for  the  generalized  poly-Weibull  distribution  with  an  arbitrary  number  of 
terms.  The  procedure  is  then  used  to  show  that  the  poly-Weibull  distribution  is 
capable  of  fitting  data  generated  from  complex  failure  processes  with  bathtub- 
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shaped  hazard  functions  better  than  the  best-fit  modified  Weibull  alternatives  in  the 
literature. 

Chapter  V  presents  the  testing  methodology  and  associated  reliability  growth 
projection  model  used  to  estimate  the  likelihood  of  achieving  a  reliability 
requirement  as  result  of  qualitative  accelerated  stress  testing.  It  is  demonstrated, 
using  published  accelerated  test  data,  that  the  proposed  model  can  assess  the 
reliability  risk  associated  with  critical  systems  or  subsystems  but  with  less  than  half 
of  the  samples  and  total  test  time  required  by  current  practices. 

To  frame  the  argument  of  accelerated  reliability  growth  modeling,  Chapter  II 
presents  an  expanded  literature  review  introducing  significant  concepts  necessary 
for  this  research  that  may  be  unfamiliar  to  the  reader.  Topics  discussed  in  this 
Chapter  include  reliability  growth  modeling,  competing  risks  analysis,  highly 
accelerated  life  testing  and  the  analysis  of  step-stress  accelerated  life  test  data. 
Chapter  VI  concludes  the  dissertation,  providing  recommendations  for  further 
implementation  of  the  model  in  design  for  reliability  contexts  and  outlines  future 
areas  of  research  to  improve  the  accuracy  and  robustness  of  the  model. 
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II.  Literature  Review 


“Six  months  in  the  lab  will  easily  save  you  a  half-day  in  the  library.  ” 

-  Ron  Kerans 


Introduction 

To  frame  the  argument  of  accelerated  reliability  growth  modeling,  this  Chapter 
presents  an  expanded  literature  review  introducing  several  major  components 
pertaining  to  the  research  that  may  be  unfamiliar  to  the  reader.  The  first  section 
provides  an  in-depth  discussion  on  the  analysis  of  time  to  event  data  when 
competing  risks  are  present.  Next,  reliability  growth  modeling  is  discussed  along 
with  several  concepts  necessary  to  link  traditional  reliability  growth  modeling  and 
accelerated  testing  to  introduce  accelerated  reliability  growth  modeling.  Then,  an 
overview  of  highly  accelerated  life  testing  (HALT)  is  presented,  with  particular 
emphasis  on  how  the  test  is  conducted,  the  purpose  of  the  test  and  how  HALT 
differs  from  other  well-known  reliability  testing  techniques.  Finally,  current 
approaches  to  modeling  time  to  failure  in  step-stress  accelerated  life  testing  are 
discussed  specifically  the  cumulative  damage  model  [18]  and  the  tampered  failure 
rate  model  [19]. 
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Competing  Risks  Analysis 


Overview 

A  competing  risk  is  an  event  whose  occurrence  fundamentally  alters  or 
altogether  eliminates  the  probability  of  observing  an  event  of  interest  [20].  As  an 
example,  the  probability  of  a  woman  developing  breast  cancer  may  become  zero  if 
she  dies  from  another  risk  factor  such  as  a  heart  attack  or  a  stroke.  Alternatively, 
the  breast  cancer  risk  for  the  same  individual  may  be  altered,  positively  or 
negatively,  after  receiving  a  diagnosis  of  lung  cancer  as  the  prescribed  medical 
treatment  could  affect  the  mechanisms  by  which  both  types  of  cancer  cells  are 
created.  Under  the  general  competing  risks  assumption  a  system  or  individual  is 
considered  to  be  at  risk  of  "failure"  from  J  risks.  The  risks  may  be  mutually 
independent  or  have  some  level  of  interdependence.  David  and  Moeschberger  [21] 
note,  however,  that  where  modeling  the  dependence  among  the  risks  was  necessary, 
authors  often  grouped  the  risks  into  categories  where  independence  among  the 
categories  could  be  assumed.  Thus  the  assumption  of  mutual  independence  among 
the  risks  has  dominated  much  of  the  competing  risks  literature. 

Under  a  given  set  of  conditions,  each  risk  competes  to  be  the  cause  of  the  failure, 
thus  the  term  risk  is  reserved  for  an  event  yet  to  occur,  while  cause  describes  the 
particular  risk  from  which  the  system  actually  failed.  Furthermore,  it  is  also 
assumed  that  systems  subject  to  competing  risks  can  fail  from  only  one  risk  and  can 
fail  only  once.  It  is  tempting  to  consider  the  survival  times  for  the  remaining/  —  1 
risks  that  were  not  the  cause  of  failure  as  randomly  right-censored,  however  this 
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type  of  non-informative  censoring  scheme  is  not  the  correct  approach  in  the 
presence  of  competing  risks.  In  non-informative  censoring,  the  only  assumption 
made  about  the  eventual  occurrence  of  an  event  of  interest  is  that  there  is  some 
positive  probability  that  the  event  will  occur  after  the  censoring  event  is  observed. 
But  in  a  competing  risks  framework  once  the  system  has  failed  due  to  risk  j  —  1, ...  ,/ 
the  probability  of  observing  the  system  fail  due  to  any  of  the  remaining/  —  1  risks  is 
altered,  and  an  informative  censoring  scheme  is  required.  Putter  et  al.  [22]  analyzed 
competing  risks  data  under  both  informative  and  non-informative  censoring 
mechanisms,  and  showed  that  a  non-informative  analysis  overestimates  the  true 
failure  probability.  Siannis  et  al.  [23]  developed  a  sensitivity  model  to  measure  the 
dependence  between  the  lifetime  of  an  individual  and  the  censoring  mechanism. 
Their  results  showed  that  the  bias  introduced  by  a  small  degree  of  dependence 
between  the  risks  can  have  a  noticeable  effect  on  the  analysis. 

Statistical  Notions  and  Notation 

In  the  competing  risks  literature,  two  statistical  notions,  represented  by  David 
and  Moeschberger  [21]  and  Crowder  [24]  dominate  the  analytical  methodology. 
Under  the  David  and  Moeschberger  notion  an  increasing  sequence  of  latent  failure 
times  is  envisioned  for  each  risk  j  =  1, ... ,/  assuming  that  j  is  the  only  risk  to  which 
the  system  is  exposed  (Figure.  2). 
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Figure  2  -  Depiction  of  competing  failure  modes  and  occurrence  times  using 

the  David  and  Moeschberger  notion 


The  notation  of  competing  risks  data  according  to  David  and  Moeschberger  is 
represented  by  the  random  variables,  Cj  and  Z,  where  Cj  is  an  indicator  representing 
the  risk  which  caused  the  failure  and  Z  =  min(7’1,  T2, ... ,  7) ),  where  7)  is  the  system’s 
theoretical  lifetime  assuming  j  is  the  only  risk  present.  If  Z  >  t,  then 

P{Z  >  t}  =  Pfa  >  t,T2>  t, ... ,  Tj  >  t)  =  Rz(t )  =  1  -  Fz(t)  (1) 


is  the  system  reliability  function.  When  mutual  independence  among  the  risks  can 
be  assumed,  (1)  becomes 

J 

Rz(t)  =  Y\Rj<t)  (2) 

i= 1 

where  Rj (t)  is  the  risk-specific  reliability  function  for  risk _/  =  1, ... Likewise,  the 


system  hazard  rate  function  for  mutually  independent  risks  is 

J 


hz(t)  = 


fzit) 

Rz(t) 


(3) 


t=l 


where  hj  (t)  is  known  as  the  cause-specific  hazard  rate  function  for  risk  j  —  1, ... ,/  in 
the  presence  of  competing  risks.  David  and  Moeschberger  [21]  further  derive  three 
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additional  statistics  commonly  used  in  the  field  of  demographics  analysis  that  may 
have  value  in  a  reliability  growth  context.  The  net  probaility  of  failure 


qj(ti,t2 )  =  1  -  exp 


l 


(4) 


measures  the  probability  of  failure  for  risk  j  assuming  it  is  the  only  risk  present. 
The  crude  probability  of  failure 


Q j(fh  £2)  —  I 

measures  the  probability  of  failure  for  risk  j  in  the  presence  of  other  competing 
risks.  Finally,  the  partial  crude  probability  of  failure 

rt  2 

Qjk(t i,t2)=  gjk(x)  exp 

is  the  probability  of  failure  for  risk  j  assuming  some  of  the  competing  risks  have 
been  eliminated  where  gjk(t)dt  and  r~k{t )  are  the  conditional  failure  probability 
and  hazard  rate  for  cause  j  in  the  absence  of  cause  k.  Equations  (4)  -  (6)  are  not 
limited  to  the  assumption  of  independent  risks,  although  when  independence  is 
assumed  the  equations  can  be  simplified  by  replacing  gj(t)  with  ry(t). 

An  alternative  notion  of  competing  risks  it  that  of  Crowder  [24],  who  modeled 
competing  risks  data  as  a  bivariate  distribution  of  the  time  to  failure,  Z,  and  an 
indicator  variable,  C,  representing  the  cause  of  failure.  Crowder  specifies  the  joint 
model  in  terms  of  the  sub-distribution  function  F(J,  t)  =  P(C  —  j,T  <  t)  or  sub¬ 
survivor  functions  R(j,  t)  —  P{C  —  j,T  >  t)  developed  by  Cox  [25]  where 

F(j,t)  + R(j,t)  =  qj.  (7) 


g j  00  exp 


- 1 
Jt, 


rz(t)dt 


dx 


(5) 
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The  sub-distribution  function  is  not  a  proper  distribution  function  as 

lim  F(J,  t)  —  q;  rather  than  unity.  Similarly,  the  sub-survivor  function  is  not  a 

proper  reliability  function  in  that  R(j,  t)  =£  P(T  >  t\C  —  j ).  The  proper  reliability 
function  is  instead  given  by  R(J,  t)/qj  where  the  proportion  qj  >  0,  j  E  1, ... ,] 
represents  the  marginal  probability  of  the  random  variable  C.  Thus  qj  — 

P(C  —  j)  —  F(j,  oo )  =  R(j,  0),  subject  to  the  constraint  YJj=1  qj  —  1.  It  follows  that 
for  failure  times  represented  by  the  continuous  random  variable  T  the  sub-density 
function,  f(j,  t)  =  —  dR(j,  t) /dt  and  the  respective  marginal  reliability  and  density 
functions  are  then  R(t )  =  Jjj=1R(j,  t)  and /(t)  =  Zy=i/0<  0-  Finally,  the  sub¬ 
hazard  function  is  expressed  as  h(t,j )  =  f(t,j)/R(t )  and  the  overall  system  hazard 
function  is  h(t)  =  HJj=1h(t,j). 

Equations  similar  to  (4)  -(6)  have  also  been  developed  using  the  Crowder 
notation.  Similar  to,  (/(iy,  t2),  the  net  probability  of  failure,  Crowder  represents  the 
distribution  of  failure  times  caused  by  risk  j  as 

U(t,j )  =  P(T  =  t\C  =j)  =  f(j,  t)/ qj.  (8) 

Alternatively,  the  distribution  of  risks  at  a  given  age  t  is  represented  by 

V(j,t)  =  (C  —  j\T  =  t)=  f(j,  t)//(t).  (9) 

Although  V (J,  t)  appears  very  different  from  the  crude  failure  probability  shown  in 
(6),  both  functions  serve  similar  functions.  Ma  and  Krings  [26]  provided  an 
excellent  comparison  of  the  David  and  Moeschberger  [21]  and  Crowder  [24] 
notions,  and  illustrated  the  conditions  under  which  both  are  equivalent. 
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Parametric  Competing  Risks  Analysis 

When  prior  knowledge  exists  on  either  the  risk-specific  system  lifetime 
distributions  or  the  underlying  failure  time  process,  parametric  approaches  to 
competing  risk  analysis  may  be  used.  These  techniques  can  reduce  the  analysis  to 
exercises  in  multivariate  statistics  wherein  the  parameter  values  may  be  estimated 
using  maximum  likelihood  or  Bayesian  estimation  techniques.  When  no  prior 
knowledge  exists  on  the  risk-specific  distributions  a  nonparametric  approach  is 
required.  These  nonparametric  approaches  are  discussed  in  the  next  section. 

In  reality,  the  risk-specific  distributions  are  rarely  known,  thus  much  of 
parametric  competing  risks  analysis  is  predicated  on  the  concept  of  theoretical  or 
latent  failure  times.  As  was  described  above,  the  latent  failure  times  parametric 
approach  presented  by  David  and  Moeschberger  assumes  that  there  exists  an 
unobservable  sequence  of  ordered  failure  times  represented  by  the  multivariate 

random  variable  Y  —  { Y u  ...  Ym . Yjl-Yjn}-  For  each  failure  mode  j,  the  sequence 

[Yj1  ...  Yjn]  indicates  the  failure  times  that  would  be  observed  if  j  were  the  only 
mode  in  the  system.  In  the  presence  of  competing  risks,  the  observable  quantities 
are  the  minimum  failure  time  Z t  —  min^n,  W-^i]  and  the  failure  cause 
indicator  Ct.  Summarizing  these  quantities,  the  observed  system  lifetime 


conditioned  on  knowing  the  cause  of  failure  is  Xt j 


Yt 


Yt  —  min(ly)  j.  Assuming 


independence  among  the  risks,  the  pdf  is  then 
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/(*() 


J 

[fviM]  ■  J~^[  5)0). 


(10) 


P  j>)  =  min  V)J  i=ij>£ 

Assuming  that  the  number  of  items  failed  due  to  cause  j  is  a  random  variable  whose 
value  can  be  represented  by  the  multinomial  distribution  function  /(%, ...  nfc)  = 

F^nf=i(p{>)  =  min  V)  jj  ,  the  likelihood  function  is  then 


n 


L  = 


n\ 


Ilf  n-i'. 


1=1 


IV  it  IV 

n  n»fc/)  ■  n  Rihi/) 


;=i 


1=1, tei 


(u) 


Defining  the  bracketed  term  as  Lit  equation  (12)  simplifies  to 


nfn(!l=l 


(12) 


Equations  (11)  and  (12)  show  that  if  the  latent  failure  times  for  each  risk  follow  a 
different  distribution,  parameter  estimation  is  accomplished  by  maximizing  each  Lj 
term  individually. 


Identifiability  Paradox 

A  disturbing  complication  exists  in  the  latent  failure  times  approach  leading  to  a 
source  of  ongoing  controversy  in  competing  risks  analysis  known  as  the 
identifiability  paradox.  Cox  [27]  first  noted  the  flaw  in  the  parametric  competing 
risks  approach  while  discussing  various  models  of  failure  times  with  two  dependent 
risks.  Tsiatis  [28]  later  generalized  the  issue  to  p  dependent  risks,  leading  Crowder 
[24]  to  rename  the  issue  as  the  Cox-Tsiatis  Impasse  -  the  issue  is  described  as 
follows.  As  discussed  in  the  beginning  of  this  section,  parametric  competing  risks 
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assumes  a  distribution  of  the  latent  failure  times  [Y1j,  Y2j  Ynj)  for  each  failure 
mode  j  under  the  assumption  that  it  is  the  only  mode  to  which  the  system  is 
subjected.  However,  since  only  the  minimum  failure  times,  Z  =  min^!,  Y\2>  -..  Y\ j\ 
are  actually  observed,  the  true  distributions  of  the  latent  failure  times  are  unknown. 
Further,  the  distribution  of  observed  lifetimes  Zit  i  —  1, 2, ... ,  /  is  completely 
determined  by  the  joint  distribution  of  the  latent  failure  times  Yjy,  however,  the 
inverse  is  not  necessarily  true.  Tsiatis  [28]  showed  that  for  any  joint  reliability 
function  with  arbitrary  dependence  between  the  risks  there  exists  another  joint 
reliability  function  with  independent  risks  that  produces  the  exact  same  sub-density 
function  f(j,  t).  Therefore  it  is  impossible  to  determine  which  model  is  correct  as 
both  will  fit  the  data  equally  well.  If  the  risks  are  independent  no  issue  exists, 
otherwise,  the  statistical  results  may  mislead  the  entire  analysis. 

Consider  an  example  in  which  two  dependent  random  variables  Y1  and  Y2  exist 
and  represent  the  failure  times  for  failure  modes  1  and  2,  respectively.  Upon 
observing  an  i.i.d.  sample  (Z,  C))  ^  Tsiatis  [28]  proved  that  for  any  distribution 

of  (Z,  C),  there  exist  independent  random  variables  Q1  and  Q2  that  provide  the  same 
distribution.  This  impasse  has  led  subsequent  research  to  focus  on  observable 
quantities  rather  than  the  joint  distributions,  thereby  estimating  the  specific  or 
"crude”  hazard  rate  l;  (t)  =  jim i  ^P(Z;-  G  [t,  t  +  At],  C;  =  l| Y  >  t),  rather  than  the 

overall  or  latent  hazard  rate  Aj  (t)  =  jim  ^  P(Zy  G  [t,t  +  At]  |Z;-  >  t). 
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Nonparametric  Competing  Risks 

Caplan  et  al.  [29]  broadly  classifies  failure  events  as  either  true  or  cause-specific 
failures  with  each  classification  having  its  own  statistical  methods  that  are  not 
necessarily  valid  for  the  other.  For  true  failures,  the  statistical  methods  are  based 
on  the  fundamental  assumption  that  all  of  the  units  under  test  would  ultimately  fail 
due  to  the  failure  mode  of  interest  were  the  study  allowed  to  continue  for  a 
sufficient  amount  of  time.  Methods  developed  to  analyze  such  data  assume  the 
underlying  cause  for  censoring  observations  is  independent  of  the  mechanism  for 
the  events  occurrence.  In  other  words,  the  survival  time  of  an  individual  (or  the 
time  at  which  a  subject  experiences  an  event)  is  assumed  to  be  independent  of  the 
mechanism  that  would  cause  the  study  to  be  censored.  In  theory,  individuals  for 
whom  the  observations  are  censored  have  an  equal  risk  of  event  of  interest 
compared  to  those  still  under  study  but  have  not  yet  observed  the  event.  This  is 
commonly  referred  to  as  non-informative  censoring  and  is  the  basis  of 
differentiation  between  the  two  commonly  used  non-parametric  estimation 
techniques  in  competing  risks  analysis:  the  complement  to  the  Kaplan-Meier 
estimator  (1  —  KM)  and  the  Cumulative  Incidence  Function  (CIF). 

Nonparametric  time-to-event  curves  are  routinely  presented  in  the  literature 
with  the  Kaplan-Meier  (KM)  product  limit  estimator  [30]  being  the  most  widely 
used.  But  this  approach  may  not  be  appropriate  when  the  analysis  is  focused  on 
time-to-first  event  in  scenarios  where  there  are  competing  events.  When  there  are 
no  competing  risks,  the  CIF  and  the  1  —  KM  estimators  produce  the  same  result.  In 
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situations  with  competing  risks,  the  cumulative  incidence  function  is  more 
appropriate.  Cheng  et  al.  [31]  notes  that  in  the  literature  the  CIF  has  been  referred 
to  as  the  cause-specific  risk,  the  crude  incidence  curve  and  the  cause-specific  failure 
probability.  Additional  references  to  CIF  described  below  are  consistent  with  the 
notation  of  Kalbfleisch  and  Prentice  [32]. 

Gooley  et  al.  [20]  and  Putter  et  al.  [22]  show  that  utilizing  the  1  —  KM 
estimator  in  the  presence  of  competing  risks  tends  to  overestimate  occurrence  rate 
of  each  event.  The  discussion  below  compares  the  two  methods,  presupposing  that 
a  population  of  items  is  at  risk  of  failure  from  two  distinct  failure  modes,  denoted 
below  as  mode  1  and  mode  2.  Further,  the  interest  is  in  estimating  the  probability 
of  failure  due  to  failure  mode  1,  vice  mode  2.  Thus,  competing  risks  are  present,  i.e. 
the  hazard  rate  functions  for  both  modes  of  failure  exist  and  the  number  of  failures 
from  the  competing  risk  will  influence  both  the  number  of  failures  and  the 
probability  of  failure  due  to  the  mode  of  interest.  The  hazard  rate  function 


h(t)  —  lim 


R(t)  -  R(t  + At) 


(13) 


At-> o  At  X  R(t) 

is  a  fundamental  relation  in  the  analysis  of  any  time-to-event  data.  At  time  t  the 
hazard  rate  is  conditioned  on  the  number  of  units  at  risk  of  failure,  thus  an 
estimator  of  h(t)  should  be  consistent  with  the  simple  ratio  of  the  number  of 
failures  divided  by  the  number  of  overall  units  under  study.  Recall  that  the  Kaplan- 
Meier  product  limit  estimator  ( KM )  is  a  nonparametric  maximum  likelihood 
estimate  of  the  reliability  function  R(t).  Conversely  the  complement  to  the  Kaplan- 
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Meier  estimator  1  —  KM  is  a  nonparametric  estimate  of  the  CDF.  Using  the  notation 
of  Kalbfleisch  and  Prentice  [32]  we  define: 
n  =  the  initial  number  of  items  at  risk 

/i  =  the  number  of  items  failed  from  the  event  of  interest  prior  to  time  tj 
rj  =  the  number  of  items  failed  from  the  competing  risk  prior  to  time  tj 
rii  =  the  number  of  items  at  risk  after  time  tj, 

and  the  total  number  of  units  at  risk  of  failure  at  any  time  t  is  expressed  as  nj  —  n  — 
(/j  +  Tj).  It  follows  that  the  respective  Kaplan-Meier  estimators  for  failure 
modes  1  and  2  are 

™lW=n(i-£) 

tj<t  t[<t 

which  shows  clearly  that  KMX  is  directly  related  to  the  hazard  rate  of  the  failure 
mode  of  interest  but  is  independent  of  the  competing  failure  mode.  Thus,  KM1(t)  is 
not  interpretable  as  the  true  probability  of  failure  mode  1  when  in  the  presence  of 
failure  mode  2. 

Alternatively,  the  cumulative  incidence  function  (CIF)  estimator  is  a  function  of 
the  hazard  rates  for  both  modes  making  CIF  the  preferred  method  of  estimating  the 
failure  probability  when  competing  risks  are  present.  Kalbfleisch  and  Prentice  [32] 
develop  a  cumulative  incidence  function  from  the  overall  Kaplan-Meier  survivor 
function  KM12(t )  =  KM1(t)  x  KM2(t)  expressed  as 

S 

CIF  if)  =  V  — f^~  X  KM12  (tj)  (15) 

Z_i  71;  —  1 
i=l 

where  s  is  the  largest  i  such  that  tj  <  t.  Gooley  et  al.  [20]  notes  that  both  1  —  KMX 
and  the  cumulative  incidence  estimator  in  (15)  are  marginal  estimates  of  the 
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probability  of  failure.  Thus,  both  estimators  will  "jump"  by  an  amount  equal  to  rT1 
after  each  occurrence  of  failure  mode  1,  until  the  first  occurrence  of  failure  mode  2. 
When  failure  mode  1  follows  failure  mode  2 

Ai -KM1  —  A CIF  X  (Xi/ni)-  (16) 

For  each  subsequent  occurrence  of  failure  mode  2  the  difference  between  the 
estimators  continues  to  expand  by  rj/n.j. 

Reliability  Growth  with  Competing  Risks 

Competing  risks  occur  frequently  in  survivability  and  reliability  analysis  and  a 
number  of  methods  have  been  proposed  for  the  analysis  of  this  type  of  data. 
However,  Ma  and  Krings  [26]  note  that  the  interaction  between  competing  risks 
analysis  and  reliability  gradually  withered  during  the  period  when  significant 
advances  were  made  in  competing  risks  analysis.  Consequently,  the  application  of 
competing  risks  analysis  in  engineering  reliability  has  fallen  behind  the  theory  of 
competing  risks  analysis. 

Corcoran  et  al.  [33]  developed  the  first  reliability  growth  projection  model  for 
estimating  reliability  after  implementing  corrective  action  in  the  final  stage  of 
development  of  an  "expensive  item."  The  Corcoran  [33]  reliability  projection  model 
is  suitable  for  use  in  cases  where  corrective  actions  are  installed  at  the  conclusion  of 
a  single  test  phase  consisting  of  N  independent  trials  where  the  number  of  trial 
outcomes  of  interest  is  a  multinomial  random  variable  with  parameters  N  (total 
number  of  trials),  p0  (unknown  initial  reliability),  and  qt  (unknown  failure 
probability  for  failure  mode  i  =  1, ... ,  k).  Since  a  multinomial  model  is  used,  the 
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equality  p0  +  Ylt=i  Pi  =  1  must  be  satisfied,  which  imposes  the  restriction  that  at 
most  one  failure  mode  can  occur  per  trial.  In  addition  to  deriving  an  exact 
expression  for  system  reliability  under  the  conditions  above,  seven  different 
estimators  are  developed  and  evaluated  for  point  estimation.  These  estimators 
were  studied  for  bias,  consistency,  conservatism,  etc.  and  ultimately  it  was  shown 
that  an  unbiased  estimate  of  the  corrected  system  could  not  be  obtained.  The 
authors  were  the  first  to  advance  the  idea  of  reducing  initial  failure  probabilities  by 
a  fractional  amount  with  consideration  to  fix  effectiveness.  Under  their  model,  the 
expected  reliability  (under  competing  risks)  at  the  end  of  the  current  test  phase  is 
given  by 

k  k 

R(N\qt )  =i?/+  'YJpiqi--'YJpi[  1  -  (1  ~qdN]  (17) 

i  i 

where  N  is  the  total  number  of  failures,  qt  is  the  failure  probability  of  failure  mode  i, 
R,  is  the  initial  reliability,  and  pt  is  the  fix  effectiveness  factor  (FEF)  for  failure  mode 
i.  Dahiya  [34]  showed  that  six  of  the  seven  estimators  initially  considered  by 
Corcoran  et  al.  [33]  possess  the  same  limiting  normal  distribution  allowing  direct 
confidence  interval  and  goodness  of  fit  procedures  for  large  samples.  Olsen  [35] 
showed  how  some  of  the  estimators  could  be  utilized  under  a  multi-stage  test 
program  and  developed  a  suitable  variant  of  the  Corcoran  model. 
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Reliability  Growth  Modeling 

Non-Accelerated  Reliability  Growth 

Introducing  a  complex  product  often  requires  a  lengthy  development  process 
during  which  it  is  expected  that  the  reliability  will  steadily  improve  based  on  testing, 
failure  mode  discovery,  root-cause  analysis  and  design  changes  or  component 
substitutions.  Testing  may  be  composed  of  many  different  types  of  tests,  each  with 
its  own  objectives.  Developmental  tests  identify  the  technical  capabilities  and 
limitations  of  proposed  designs  and  ensure  sufficient  design  maturity  is  achieved 
prior  to  operational  testing.  In  operational  testing  the  focus  is  on  demonstrating 
that  the  design  is  suitable  for  its  intended  use  in  a  realistic  operational  environment. 
Entering  operational  test  with  an  immature  design  often  results  in  continued 
debugging  into  the  early  life  of  the  product  after  it  has  been  released  to  the  market, 
usually  at  much  greater  cost  than  if  the  fault  were  discovered  in  developmental 
testing.  Therefore,  it  is  desirable  to  model  the  improvement  in  reliability  over  time 
to  (1)  forecast  the  length  of  the  development  process,  (2)  ensure  proper  allocation 
of  testing  resources  and  (3)  estimate  the  reliability  upon  entry  into  market. 

The  initial  product  of  a  reliability  growth  model  is  the  idealized  reliability 
growth  planning  curve  (Figure  3).  Created  early  in  the  development  process,  the 
idealized  curve  is  a  roadmap  to  baseline  the  reliability  growth  progress  within  a 
single  test  phase.  When  testing  continues  across  test  phases,  multiple  idealized 
curves  are  developed.  An  indication  that  the  reliability  is  lagging  (assessments 
below  the  idealized  curve)  may  signal  the  need  to  reallocate  test  resources  or 
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reevaluate  the  testing  strategy  to  ensure  the  reliability  achieves  the  requirement 
with  a  specified  level  of  confidence.  Many  models  plot  reliability  growth  as  an 
increase  in  mean  time  between  failures  ( MTBF ),  or  a  decrease  in  the  failure  rate 
(MTBF-1)  against  cumulative  test  time  across  all  units  at  risk  [36]. 


Planned  Growth 
Assessed  Growth 


Reliability 


Test  Phase  1  Test  Phase  2  Test  Phase  3 

Figure  3  -  Idealized  reliability  growth  planning  curves  for  three  test  phases  showing 
the  difference  in  planned  and  assessed  reliability  (Source:  MIL-HDBK-189C,  p.  11) 

When  many  modular  tests  and  debugging  efforts  are  proceeding  in  parallel, 
growth  may  be  measured  as  the  rate  at  which  defects  are  found  and  corrected  in  a 
testing  interval.  In  principle,  reliability  growth  models  may  also  allow  for  decreases 
in  reliability.  For  example,  a  product’s  reliability  may  be  adversely  affected  as  a 
result  of  substituting  a  cheaper  but  less  reliable  component.  Our  focus  here,  though, 
is  on  positive  reliability  growth  which  improves  product  reliability. 

According  to  Ebeling  [37]  the  earliest  developed  and  most  frequently  used 
reliability  growth  model  is  that  of  Duane  [38]  who  postulated  from  his  empirical 
observations  that  for  products  under  development,  plots  of  cumulative  failure  rate 
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versus  cumulative  operating  hours  on  a  log-log  scale  were  approximately  linear 
with  negative  slope.  Crow  [39]  derived  the  stochastic  basis  for  the  Duane  model, 
referred  to  in  the  literature  as  a  power  law  process  (PLP)  or  informally  as  a  Weibull 
process.  The  Crow  [39]  model  is  a  type  of  nonhomogeneous  Poisson  process 
(NHPP)  with  a  Weibull  hazard  intensity  function  h(t)  =  where  X  > 

0  and  /?  >  0  are  interpreted,  in  the  Crow  model  context,  as  the  initial  reliability  of 
the  system  and  the  growth  parameter,  respectively.  For  /?  <  1,  the  hazard  is 
decreasing  (reliability  is  increasing)  while  /?  >  1  implies  an  increasing  hazard  and 
/?  =  1  means  that  both  the  reliability  and  the  hazard  rate  are  unchanged.  Duane 
[38]  observed  /?  values  in  the  neighborhood  of  0.50  while  the  Army  Materiel 
Sustainment  Analysis  Activity  (AMSAA)  [9]  suggests  values  in  the  range  (0.30,  0.75), 
depending  on  the  level  of  commitment  to  reliability  improvement  and  the  type  of 
system  under  development.  The  rigorous  development  of  the  PLP  allows  estimation 
and  prediction  in  both  the  maximum  likelihood  estimation  (MLE)  and  Bayesian 
contexts  and  has  made  the  approach  popular  in  the  literature. 

However,  the  PLP  and  other  models  based  on  the  Duane  postulate  are  an 
idealization  of  the  true  underlying  Test-Analyze-And-Fix  (TAAF)  failure  process. 

Sen  [40]  notes  that  these  models  estimate  the  improvement  in  reliability  based  on 
the  number  of  failure  modes  discovered  without  accounting  for  fixes  or  design 
changes.  Various  reliability  growth  models  have  been  developed  since  the  Duane 
Postulate  to  plan,  track,  and  project  the  reliability  improvement  of  a  system 
throughout  the  development  process.  Planning  models  are  used  to  develop  the 
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idealized  and  planned  growth  curves  to  address  program  schedules,  required 
testing  resources  and  the  realism  of  achieving  the  reliability  requirement  with  a 
desired  level  of  confidence  within  the  allocated  test  program.  Tracking  models 
gauge  the  progress  of  the  reliability  effort  based  on  test  results.  Projection  models 
indicate  the  anticipated  reliability  at  some  future  time  based  on  achievement  to  date 
and  engineering  assessments.  Hall  [10]  provides  one  of  the  most  detailed  and 
comprehensive  reviews  of  reliability  growth  planning,  tracking  and  projection 
models  for  both  continuous  use  and  discrete  use  systems  found  anywhere  in  the 
literature.  The  review  comprises  a  synopsis  of  over  80  papers  covering  planning 
models,  tracking  and  projection  models.  Further,  numerous  reliability  growth 
surveys/handbooks  and  thirty-six  other  papers  covering  theoretical  results, 
simulation  studies,  real-world  applications,  personal  perspectives,  international 
standards,  or  related  statistical  procedures  are  also  included. 

Accelerated  Reliability  Growth 

As  opposed  to  the  reliability  growth  models  discussed  by  Hall  [41]  in  which 
growth  is  based  on  a  decreasing  rate  of  failure  mode  discovery  in  successive  test 
intervals,  step-stress  testing  results  in  an  increasing  number  of  failure  modes  as 
testing  progresses  to  higher  stresses.  Further,  step-stress  testing  is  a  continuous 
process  with  no  discernible  intervals  and  the  removal  of  each  failure  mode 
discovered  does  not  yield  the  same  growth  at  the  non-accelerated  use  stress. 
Therefore,  it  is  claimed  that  an  improved  estimate  of  the  system’s  field  reliability 
can  be  obtained  by  analyzing  HALT  failure  data  as  a  test-fix-test  process  and 
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distinguishing  between  the  modes  of  failure  that  occur  in  a  test  through  the 
incorporation  of  competing  risks  in  the  analysis. 

Upon  starting  a  traditional  reliability  growth  test,  it  is  assumed  that  few  failure 
modes  remain  in  the  system  for  which  the  root  cause  has  not  been  identified  and 
understood.  However,  even  when  this  assumption  is  met  it  is  not  uncommon  for 
growth  testing  to  require  several  months  to  attain  the  reliability  goal  with  the 
required  confidence  level.  Thus,  accelerated  test  methods  are  increasingly  utilized 
as  part  of  a  reliability  growth  strategy  to  reduce  the  overall  test  duration.  Several 
widely  used  design  improvement  techniques  [12,  42, 11, 13,  43]  expose  early 
product  designs  to  accelerated  environmental  and  repetitive  use  stresses  to  quickly 
force  latent  defects  to  manifest  themselves  thereby  ensuring  a  highly  reliable 
product  is  fielded  quickly.  The  defining  characteristic  of  these  tests  lies  in  the 
underlying  assumption  that  simply  meeting  the  validation  requirement  is 
insufficient  for  ensuring  high  field  reliability  as  the  requirement  rarely  encompasses 
the  full  stress  envelope  a  product  will  actually  encounter  over  its  service  life.  Thus, 
in  many  instances  failure  modes  exposed  by  these  techniques  are  not  representative 
of  field  experience,  but  reflect  brief  high-level  stress  states  often  not  considered. 

In  a  reliability  growth  test  it  is  assumed  that  (1)  a  set  of  latent  failure  modes 
exists  in  each  system  under  test,  (2)  the  test  profile  will  activate  a  sufficient  number 
of  these  latent  modes  during  the  allocated  test  time  and  (3)  corrective  action 
removes  a  fraction  of  each  discovered  mode’s  failure  intensity.  In  other  words,  if  the 
mode  specific  failure  rate  due  to  mode  i  is  represented  by  A*,  the  overall  system 
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failure  rate  is  Y.i  If  each  failure  mode  receives  corrective  action  that  reduces  by 
the  proportion  p,  one  would  conclude  that  fixing  any  failure  mode  reduces  the 
overall  system  failure  rate  by  1  —  p). 

Often  the  test  profile  involves  the  application  of  constant  stresses  chosen  to 
represent  field  conditions.  Thus  the  proper  selection  of  stressors  and  stress  levels  is 
critical  to  achieving  adequate  reliability  growth.  If  the  stress  profile  is  a  function  of 
time,  however,  consideration  must  be  given  to  the  stress  level  at  which  a  failure 
mode  is  discovered.  Under  progressive  stress  or  step-stress  profiles  it  is  possible  to 
discover  failure  modes  at  high  stress  levels  that  would  never  be  activated  at  the  use 
stress  level  regardless  of  the  exposure  time.  Thus,  implementing  corrective  action 
on  these  failure  modes  may  produce  little  to  no  actual  reliability  growth. 

Highly  Accelerated  Life  Testing  (HALT) 

Highly  accelerated  life  testing  originated  in  the  electronics  industry  when  Hobbs 
began  to  use  the  term  in  1988  [44].  Though  the  mindset  and  process  of  HALT  goes 
back  several  decades,  Hobbs  coined  the  term  to  describe  a  sequential  testing 
process  wherein  products  are  subjected  to  a  variety  of  stresses  to  identify  the  weak 
links  in  the  design.  More  formally,  HALT  involves  the  application  of  individual  or 
simultaneous  stressors  at  levels  elevated  beyond  those  experienced  in  either  the  use 
environment  or  in  traditional  accelerated  life  tests.  During  testing,  failure  modes 
are  discovered  and  removed  from  the  design  through  corrective  action  in  an 
iterative  process  which  expands  the  product’s  operating  margin.  Once  the  test  is 
complete,  the  product  will  have  reduced  random-failure  probabilities  and  longer 
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lifetimes  for  wear-out  in  the  designed  usage  environment.  It  should  be  noted, 
however  that  HALT  does  not  specifically  address  wear-out  failure  modes. 

According  to  McLean  [12],  HALT  proceeds  by  first  applying  individual  stressors 
that  are  increased  step-wise  at  periodic  intervals  until  a  failure  event  occurs  (i.e.  no 
censoring).  No  standard  exists  to  direct  what  stressors  should  be  included  for  a 
given  product  type.  However,  McLean  notes  that  experience  has  led  practitioners  to 
adopt  an  ordered  set  of  stressors  (Table  1)  as  standard  to  HALT  which  represents 
what  is  used  in  the  majority  of  tests  conducted  to  date.  As  the  goal  of  HALT  is  to 
remove  as  many  failure  modes  as  possible,  practitioners  are  not  limited  to  the 
stressors  listed  in  Table  1  nor  are  they  required  to  ensure  that  each  stressor  is 
included.  Product  specific  stresses  are  often  included  either  independently  or  in 
combination  with  another  test. 


Table  1  -  HALT  stress  sequence 

Stressor  Notes 


%  Total  Failure 
Modes 


Cold  Step  Stress  14% 

Hot  Step  Stress  17% 

Rapid  Thermal  Transitions  >60C°/minute  4% 

Vibration  Step  Stresses  6  DoF  random  vibrations  45% 

Rapid  Thermal/Vibration  20% 


When  a  failure  event  occurs  in  HALT,  a  key  distinction  is  made  in  identifying  the 
event  as  either  a  hard  failure,  an  operating  limit  failure  or  a  destruct  limit  failure. 
This  distinction  guides  the  HALT  process,  indicating  whether  testing  should  proceed 
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at  higher  levels  of  stress  or  if  sufficient  progress  has  been  made  and  the  test  should 
be  stopped.  Hard  failures  result  from  damage  to  the  system’s  physical  architecture 
requiring  corrective  action  to  restore  performance  to  specification  levels.  Operating 
limit  (OL)  failures  occur  when  one  or  more  performance  specifications  aren’t  being 
met  and  cannot  be  brought  back  to  within  specification  levels  without  reducing  the 
stress  level.  It  should  be  noted  that  the  occurrence  of  an  OL  failure  does  not 
necessarily  imply  the  occurrence  of  a  hard  failure.  Destruct  limit  (DL)  failures,  also 
known  as  the  fundamental  limit  of  the  technology  (FLT),  describe  the  occurrence  of 
unavoidable  failure  modes  due  to  the  physical,  chemical  or  structural  limits  of  an 
item.  Removing  these  failure  modes  typically  requires  a  material  substitution  or 
significant  structural  change  that  is  either  physically  impossible  or  financially 
infeasible.  Examples  of  DL  failures  include  the  melting  or  phase  change  of  a 
material,  catastrophic  mechanical  failure  due  to  vibrational  stresses  and  dielectric 
breakdown  due  to  voltage  overstress.  DL  failures  modes  would  not  be  considered 
representative  of  field  experience  and  are  useful  only  for  driving  system 
performance  to  the  highest  levels  possible. 

A  product’s  operating  margin  is  the  stress  range  in  which  the  system  can  operate 
for  a  specified  mission  duration  without  experiencing  an  operational  limit  failure. 
Figure  4  illustrates  the  relationship  between  a  product’s  specification  and  its 
operating  and  destruct  limits  before  and  after  a  HALT.  Comparing  the  top  and 
bottom  graphics  in  Figure  4  illustrates  clearly  the  outward  shift  in  the  product’s 
operating  margin  after  the  HALT.  Additionally,  Figure  4  suggests  that  the  upper  and 
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lower  operating  and  destruct  limits  of  a  system  are  random  quantities  whose  values 
can  be  represented  by  normal  distributions  with  respective  mean  fi  and  variance  o. 

maxSj 

These  parameters  could  then  be  estimated  by  the  sample  means  fii  —  — — - and 

X^lmaxSj-g] 

sample  variances  0)  —  — — -  - of  the  resulting  maximum  operating  limit  stress 

levels  attained  in  testing  each  of  the  samples. 
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Figure  4  -  Illustration  of  product  specification,  operating  and  destruct  limits  before 

and  after  highly  accelerated  life  testing 

Despite  superficial  similarities  in  terminology,  HALT  is  disjoint  from  accelerated 
life  testing  (ALT)  in  several  key  areas.  Whereas  physical  models  have  been 
developed  [45]  to  compute  acceleration  factors  that  allow  the  inclusion  of  ALT  test 
results  into  a  reliability  calculation,  no  such  models  exist  for  the  inclusion  of  HALT 
test  data  due  (1)  to  their  inability  to  generalize  across  larger  stress  extrapolations, 
(2)  the  potential  for  complex  stress  profiles  and  (3)  changes  in  the  product’s  design 
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configuration.  Nevertheless,  practitioners  have  sought  ways  to  predict  field 
reliability  from  the  results  of  HALT.  Silverman  [46]  discussed  the  difficulties  in 
analyzing  HALT  data  using  physics  of  failure  models  and  acceleration  factor 
calculations  due  primarily  to  limited  data  sharing  across  product  categories 
allowing  for  more  generalized  approaches.  Meeker  [17]  expressed  concerns  with 
the  complexity  of  HALT  stress  profiles,  advising  practitioners  to  avoid  analyzing 
HALT  data  using  techniques  derived  from  ALT.  However,  Meeker’s  discussion 
targeted  tests  in  which  multiple  stresses  were  applied  simultaneously  and  did  not 
address  simpler  single  stress  test  designs  which  are  more  common. 

Sequential  test  procedures  have  been  developed  to  utilize  the  results  of  HALT  as 
a  guide  for  follow-on  tests  to  reduce  the  error  involved  in  extrapolations  across 
stress  levels.  Edson  [47]  introduced  Calibrated  Accelerated  Life  Testing  (CALT)  as  a 
General  Motors  corporate  standard  to  evaluate  the  reliability  of  products  subject  to 
HALT  by  adding  test  samples  at  key  stress  levels.  The  choice  of  the  number  of 
additional  samples  and  stress  levels  to  be  tested  in  CALT  are  based  on  the  maximum 
stress  derived  from  the  HALT.  Bhote  and  Bhote  [11]  developed  a  nine-step  life 
prediction  modeling  process  known  as  multi-environment  over-stress  testing 
(MEOST).  This  method  bypasses  the  independent  stress  tests  and  applies  all  of  the 
stresses  at  once.  Reliability  estimation  with  MEOST  is  accomplished  after  observing 
a  full  year  of  service  life  to  produce  a  useful  estimate  of  the  product’s  reliability 

Brand  and  McLean  were  awarded  three  patents  [48,  49,  50]  comprising  the 
proprietary  HALT  AFR  Calculator®  (AFR)  model  currently  managed  by  the  testing 
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firm  Ops  A  La  Carte®.  The  HALT  AFR  Calculator®  is  the  most  prominent  HALT 
analysis  model  to-date  utilizing  linear,  quadratic  and  exponential  time  acceleration 
techniques  to  estimate  a  product’s  relative  life  expressed  as 

MTBF2 


R  = 


(18) 


MTBF1 

where  MTBF1  is  the  field  use  reliability  for  the  original  unit  subjected  to  an  initial 
HALT  and  MTBF2  is  the  field  reliability  for  the  unit  subjected  to  a  second  HALT  after 
implementing  corrective  action  to  remove  the  failure  discovered  during  the  first 
HALT.  To  utilize  the  HALT  AFR  Calculator  it  is  assumed  that  1)  an  estimate  of  the 
system’s  initial  field  reliability  is  known  prior  to  testing,  2)  that  two  rounds  of  HALT 
comprised  of  k  independent  stress  tests  have  been  performed  and  3)  the  stress 
limits  for  both  rounds  of  HALT  are  identical.  Applying  only  the  individual  tests 
listed  in  Table  1,  the  AFR  estimates  relative  life  as  a  function  the  ratio  of  the  times  to 


first  failure  in  HALT1  and  HALT2  as 


R  = 


f(tA  i) 


(19) 


/  (L42) 

where  tA1  and  tA2  are  assumed  to  be  realizations  of  independent  exponential 
distributions  and  /(  )  implies  either  a  linear,  quadratic  or  exponential  function  of 
the  failure  times,  depending  on  which  of  the  three  patents  is  considered.  Applying 
an  appropriate  acceleration  factor  (i.e.  Arrhenius,  Inverse  Power  Law,  Eyring,  etc.) 
gives  the  ratio  of  times  to  first  failure  in  the  un-accelerated  time  scale 

/  (Fai)  f  (tfi) 


R  —  AF 


f  (L42)  f^Fl) 


(20) 
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Taking  In  [r]  gives  R*  —  /( tF  —  tF  )  =  /(AtF),  for  each  of  the  k  individual  stressor 
tests,  an  independent  R*k  is  obtained.  The  point  estimate  for  the  relative  life  R  is 
then  found  from  the  simple  average  of  these  values  and  confidence  limits  are 
derived  assuming  that  each  Rk  is  approximately  normally  distributed. 

Though  McLean  [51]  has  presented  the  results  of  the  AFR  Calculator  to  the 
accelerated  test  and  reliability  community,  the  manner  in  which  the  linear, 
quadratic  and  exponential  acceleration  techniques  are  utilized  are  unpublished. 

Step  Stress  Accelerated  Life  Testing  Data  Analysis 

Step-stress  accelerated  life  testing  (SSALT)  involves  testing  items  under 
successively  higher  levels  of  stress  to  drive  failures  quickly  and  reduce  the  overall 
testing  time.  The  SSALT  strategy  avoids  the  problem  inherent  to  constant  stress 
accelerated  life  testing  of  selecting  the  stress  levels  to  generate  a  sufficient  number 
of  failures  to  enable  statistical  inference.  In  SSALT,  an  initial  load  is  applied  at  time 
t  —  0,  for  a  specified  duration,  after  which  the  load  is  increased  step-wise  as  time 
thresholds  are  reached.  Testing  stops  when  (1)  a  sufficient  number  of  items  on  test 
have  failed,  (2)  a  limit  stress  has  been  reached  or  (3)  the  total  time  available  for 
testing  has  been  exhausted.  The  time  duration  that  the  units  dwell  at  a  given  stress 
level  may  be  based  on  the  time  required  to  conduct  system  diagnostics  or  on  some 
optimality  criterion.  Bai  [52]  derived  the  optimal  step  time  based  on  minimizing  the 
asymptotic  variance  of  the  maximum-likelihood  estimator  of  the  mean  life  at  a 
design  stress.  Laio  [53]  used  the  asymptotic  variance  of  the  estimated  lOOpth 
percentile  of  the  product's  lifetime  distribution.  In  a  HALT,  equivalent  time 
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increments  are  used  and  the  length  of  each  interval  is  based  on  the  time  required  to 
conduct  the  necessary  system  performance  diagnostics  at  each  stress. 

Estimating  the  lifetime  of  products  at  a  use  stress  is  substantially  more  difficult 
in  products  exposed  to  SSALT  as  the  cumulative  exposure  at  each  stress  level  must 
be  accounted  for  simultaneously.  Two  models  commonly  used  for  this  purpose  are 
the  cumulative  hazards  model  by  Nelson  [54]  and  the  tampered  failure  rate  model 
by  Bhattacharyya  and  Soejoeti  [19].  Some  variations  and  generalizations  of  these 
models  are  presented  by  Kececioglu  [55]  and  Zhao  et  al.  [56].  In  this  literature 
review,  only  the  Cumulative  hazards  model  is  discussed.  It  should  be  noted, 
however,  that  when  the  baseline  failure  distribution  is  exponential,  the  results 
obtained  by  both  models  coincide. 

Cumulative  Hazards  Model 

Consider  the  ordered  sequence  of  increasing  stresses  (i  =  1, 2, ... )  where  the 
ratio  si+1/sj  is  constant  with  respect  to  time  and  the  product’s  underlying  life 
distribution  at  use  stress  s0  is  Weibull.  If  an  accelerated  stress  sx  is  applied,  the 
fraction  of  units  failing  up  to  time  t  is  given  by 

/qCt;  Si)  =  1  -  exp{-( Ks^tY)  .  (21) 

where  the  Weibull  scale  parameter  a  —  Ks^n  reflects  an  inverse  power  law  life- 
stress  relationship  with  parameters  K  and  n.  If  the  stress  level  is  changed  to  s2  (21) 
is  no  longer  appropriate  to  model  system  lifetime  as  it  does  not  incorporate  the 
cumulative  exposure  at  both  stress  levels.  Thus,  to  analyze  the  data  from  a  step- 
stress  test  perspective,  a  cumulative  exposure  model  is  needed. 
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The  cumulative  hazards  model  [18]  relates  the  life  distribution  of  units  under 
test  at  one  stress  level  to  the  life  distribution  at  another  stress  level.  It  is  assumed 
that  the  remaining  life  of  the  test  units  depends  only  on  the  cumulative  exposure  the 
units  have  experienced  up  to  the  current  time  and  is  independent  of  how  the 
exposure  was  accumulated.  Moreover,  since  the  units  are  held  at  a  constant  stress 
at  each  step,  the  surviving  units  will  fail  according  to  the  distribution  at  the  current 
step,  but  with  a  starting  age  corresponding  to  the  total  time  accumulated  prior  to 
the  beginning  the  current  step.  Thus,  the  probability  of  units  failing  during  the  time 
interval  (0,  tx)  under  stress  sx  not  having  experienced  any  other  stresses  may  be 
described  by  (21).  After  time  the  probability  that  the  surviving  units  will  fail 
during  (t1;  t2),  while  exposed  to  stress  level  s2  is  equivalent  to  the  probability  that 
the  units  would  fail  while  at  s2  after  accumulating  t  —  tx  plus  an  equivalent  age  ex  to 
account  for  the  exposure  at  sx  expressed  as 

^2(t;s2)  =  1  -  exp  [-  (Ks?((t  -  tj)  +  (22) 

The  equivalent  age  ex  is  the  time  at  which  the  CDF  for  s2  is  equal  to  the  CDF  for 
sx  after  an  exposure  of  tx,  or 
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This  process  may  be  used  recursively  to  for  an  arbitrary  number  of  stress  levels 
where,  in  general  the  equivalent  age  is  expressed  as 


Conclusion 

This  literature  review  has  provided  a  framework  for  merging  accelerated 
testing  techniques  with  reliability  growth  methods  by  introducing  several  major 
components  pertaining  to  the  research  that  may  have  been  unfamiliar  to  the  reader. 
The  Chapter  has  provided  a  synopsis  of  the  research  accomplished  in  the  fields  of 
reliability  growth  for  complex  systems,  competing  risks  analysis  and  accelerated  life 
testing  data  analysis.  The  parametric  and  nonparametric  analysis  of  time  to  event 
data  when  competing  risks  are  present  was  first  discussed.  Two  statistical  notions 
were  presented  as  was  the  non-identifiability  problems  that  can  result  when  the 
risks  cannot  be  assumed  to  be  mutually  independent.  Reliability  growth  modeling 
was  then  discussed  as  were  several  concepts  necessary  to  link  traditional  reliability 
growth  modeling  and  accelerated  testing  to  introduce  accelerated  reliability  growth 
modeling.  Finally,  the  cumulative  damage  model  [18]  was  discussed  as  an  approach 
to  model  time  to  failure  in  step-stress  accelerated  life  testing.  In  the  following 
Chapters,  much  of  the  information  presented  in  this  literature  review  is  utilized  to 
develop  the  reliability  growth  projection  model  to  translate  failure  data  obtained 
from  a  qualitative  accelerated  life  test  into  to  an  estimate  of  system  reliability  after 
implementing  corrective  action. 
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III.  Bridging  the  Gap  between  Quantitative  and  Qualitative  Accelerated  Life 

Tests 
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Abstract  -  Test  planners  have  long  sought  the  ability  to  incorporate  the  results  of  highly 
accelerated  life  testing  (HALT]  into  an  early  estimate  of  system  reliability.  While  case-studies 
attest  to  the  effectiveness  of  HALT  in  producing  reliable  products,  the  capability  to  translate 
the  test's  limited  failure  data  into  a  meaningful  measure  of  reliability  improvement  remains 
elusive.  Further,  a  review  of  quality  and  reliability  literature  indicates  that  confusion  exists 
over  what  defines  a  highly  accelerated  life  test  and  how  HALT  differs  from  quantitative 
accelerated  life  testing  [QALT]  methods.  Despite  many  authors  making  a  clear  distinction 
between  qualitative  and  quantitative  accelerated  life  tests,  an  explanation  as  to  why  this 
delineation  exists  cannot  be  found.  In  this  paper,  we  consider  an  exemplary  HALT  composed 
of  a  single  stressor  to  show  that  the  HALT  philosophy  precludes  the  estimation  of  a  system’s 
hazard  rate  function  parameters  due  to  the  test's  fix  implementation  strategy.  Four  common 
accelerated  failure  data  analysis  methods  are  highlighted  to  show  their  limitations  with 
respect  to  estimating  reliability  from  HALT  data.  Finally,  a  way  forward  for  future  research 
is  provided  for  improving  the  parameter  estimation  in  follow-on  testing. 

Index  Terms  -  competing  risks,  highly  accelerated  life  testing,  reliability  growth,  step- 
stress  accelerated  life  testing 


Introduction 

Test  planners  have  long  sought  the  ability  to  incorporate  the  results  of  highly 
accelerated  life  testing  (HALT)  into  an  early  estimate  of  system  reliability.  While 
case-studies  [15, 16, 14]  attest  to  the  effectiveness  of  HALT  in  producing  reliable 
products,  the  capability  to  translate  the  test’s  limited  failure  data  into  a  meaningful 
measure  of  reliability  improvement  remains  elusive.  Further,  a  review  of  quality 
and  reliability  literature  indicates  that  confusion  exists  over  what  defines  a  highly 
accelerated  life  test  and  how  HALT  differs  from  quantitative  accelerated  life  test 
(QALT)  methods.  Despite  authors  making  distinctions  between  qualitative  and 
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quantitative  accelerated  life  tests,  a  clear  explanation  as  to  why  this  delineation 
exists  cannot  be  found  in  the  literature. 

The  current  paper  introduces  the  purpose  and  process  of  a  HALT  before 
presenting  an  exemplary  single-stressor  highly  accelerated  life  test  to  discuss  the 
test’s  data  structure.  Then,  approaches  for  characterizing  product  reliability  are 
examined  to  highlight  their  limitations  with  respect  to  parameter  estimation  from 
HALT  data.  Finally,  avenues  for  future  research  to  utilize  HALT  data  in  follow-on 
testing  are  then  suggested. 

Highly  Accelerated  Life  Testing 

HALT  originated  in  the  electronics  industry  when  Hobbs  [44]  began  to  use  the 
term  to  describe  a  sequential  testing  process  capable  of  quickly  identifying  the  weak 
links  in  a  product’s  design.  More  formally,  HALT  is  an  iterative  test-fix-test  [57] 
reliability  growth  process  (Figure  5)  wherein  failure  modes  are  discovered  and 
removed  as  result  of  applying  one  or  more  stressors  beyond  the  levels  experienced 
in  either  the  use  environment  or  traditional  accelerated  life  tests  [12].  By  subjecting 
products  to  such  extreme  environments  many  hard-to-find  failure  modes  that 
would  otherwise  go  undetected  in  traditional  reliability  tests,  can  be  discovered  and 
removed  quickly.  The  HALT  methodology  is  motivated  by  the  tendency  of  nearly  all 
products  to  experience  one  or  more  excursions  outside  of  the  expected  usage 
environment  over  its  lifetime.  Such  excursions  can  occur  if  the  total  stress  exposure 
is  not  anticipated  in  product  design  (e.g.  transportation  stresses)  or  if  the  product  is 
misused  after  fielding. 
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Increase  Stress 


Induce  Failure 


Corrective  Action 


Failure  Analysis 


Figure  5  -  Four-step  highly  accelerated  life  testing  process 


HALT  improves  a  product’s  reliability  by  expanding  its  operating  margin  (Figure  6) 
as  failure  modes  are  discovered  and  removed,  creating  a  "safety  net"  to  withstand 
these  excursions. 

Lower  Lower  Upper  Upper 

Destruct  Opcr.  p  ,  Opcr.  Destruct 


Prior  to  HALT 

Lower  Lower  Upper  Upper 


After  HAM 

Figure  6  -  Relationship  between  product  operating  and  destruct  limits 


No  formal  specification  exists  to  describe  how  HALT  should  be  conducted,  but  a 
procedure  described  by  Mclean  [12]  and  used  in  industry  is  considered  common 
practice  [58,  59].  In  this  procedure,  four  to  eight  prototypes  are  exposed  to  an 
ordered  stress  regimen  (Table  2).  However,  since  the  goal  of  HALT  is  to  discover  as 
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many  latent  failure  modes  as  possible,  test  plans  are  neither  limited  to  these 
stressors  nor  are  they  required  to  ensure  that  each  stressor  is  included.  Often 
product  specific  stresses  (e.g.  voltage/frequency  margining)  are  added  to  or 
combined  with  these  "standard”  tests. 


Table  2  -  McLean  HALT  stress  regimen 


Stress  Order 

Stressor  Type 

Application  Notes 

1 

Cold  Step  Stress 

2 

Hot  Step  Stress 

3 

Rapid  Thermal 
Transitions 

>  60C°/minute  rate  of 
temperature  change 

4 

Vibration  Step  Stresses 

Simultaneous  six  degree  of 
freedom  random  vibrations 

5 

Combined  Rapid  Thermal 
Transitions/Vibration 

Vibration  applied  at  thermal 
limits  identified  in  Steps  1  and  2 

During  HALT,  prototypes  are  electronically  monitored  to  detect  any  degradation 
in  performance  or  loss  of  function.  When  a  failure  occurs  or  a  performance  measure 
cannot  be  brought  back  within  specification  limits  via  stress  reduction,  corrective 
action  is  required  to  alter  the  system’s  architecture  and  remove  the  design 
weakness.  As  a  result,  the  improved  system  can  withstand  higher  stress  levels  prior 
to  failure.  This  test-fix-test  process  continues  iteratively  as  newer  designs  are 
stressed  up  to  the  physical  or  chemical  limits  of  the  unit  (e.g.  melting,  dielectric 
breakdown).  At  this  fundamental  limit,  corrective  action  requires  a  financially 
prohibitive  design  change,  thereby  concluding  the  test  and  finalizing  the  design. 
Upon  release,  the  updated  product  will  have  reduced  infant  mortality  and  random 
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failure  probabilities  along  with  longer  lifetimes  with  respect  to  wear-out  in  the 
designed  usage  environment  as  shown  in  Figure  7. 

—  Before  HALT 


Infant  Useful  Wear 

Mortality  Life  Out 


Figure  7  -  Effect  of  highly  accelerated  life  testing  on  bathtub  curve 

Example  HALT  Scenario  and  Data  Structure 

Meeker  et  al.  [17]  discuss  the  potential  pitfalls  of  interpreting  failure  data 
produced  from  accelerated  tests  with  complex  stress  profiles  and  multiple  stressors. 
Here,  an  exemplary  HALT  composed  of  a  single  stressor  is  posited  to  illustrate  that 
the  HALT  process,  not  the  stress  profile,  precludes  estimating  a  system’s  hazard  rate 
function  parameters.  Consider  a  stress  profile  where  HALT  testing  begins  with  the 
application  of  an  initial  stress  level  s0  and  increases  step-wise  (Figure  8)  according 
to  the  piece-wise  right-continuous  function 

s(t)  =  s0  +  8  ■  Z(t )  (25) 

where  Z(t)  =  max{i\ <  t,  i  —  0, 1, ...  ,  M},  8  is  the  fixed  step-up  stress  increment, 
and  [t j,  ri+1)  represents  the  time  interval  in  which  the  ith  stress  level,  Sj,  is  applied. 
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The  sample  set  of  units  at  risk  in  this  exemplary  HALT  is  comprised  of  m  identical 
prototypes,  sharing  an  initial  design  configuration  Cj  —  c0.  Like  most  developed 
products,  the  prototypes  are  at  risk  of  failure  from  multiple  causes  where  each 
cause  is  itself  the  result  of  a  flaw  in  the  design  activating  in  response  to  a  given 
input.  Some  of  these  flaws  may  be  found  easily  in  that  the  range  of  stress  inputs 
which  are  likely  to  expose  them  is  large.  However,  many  more  flaws  will  be 
relatively  difficult  to  find  since  only  a  small  number  of  stress  inputs  will  reveal  their 
presence. 


s,  -- 
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aj 

> 

aj 
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Time 


Figure  8  -  Example  step-stress  profile 


During  [Tj,  ri+1),  i  —  1, 2  ... ,  M  an  unknown  number  of  latent  failure  modes,  denoted 
by  N,  compete  to  be  the  cause  of  failure.  Not  all  of  the  N  overall  modes  are  likely  to 
be  discovered  as  result  of  exposure  to  a  single  stressor  (Figure  9),  thus  for  the 
single-stressor  HALT  considered  here  we  define  J  <  N  as  the  number  of  failure 
modes  ultimately  discovered  during  the  test. 
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Figure  9  -  Fraction  of  failures  discovered  by  common  HALT  stressors  (McLean,  2008) 

This  framework  is  akin  to  a  serial  arrangement  of  the  failure  modes  where  the 
occurrence  of  mode  j  E  1, ...  ,J  at  stress  level  st  causes  overall  system  failure  and 
interrupts  the  test  at  clock  time  t;  where  Tt  <  tj  <  ri+1.  Implementing  corrective 
action  results  in  a  new  design  configuration  wherein  failure  mode  j  no  longer  exists, 
but  additional  failure  modes  not  contained  in  the  previous  design,  may  have  been 
created.  Testing  resumes  on  the  updated  design  at  stress  st  once  the  corrective 
action  has  been  implemented  on  the  unit  under  test  as  well  as  any  remaining 
untested  units.  Consequently,  the  HALT  test-fix-test  strategy  allows  for  a  single 
observation  of  failure  mode  j  prior  to  testing  the  new  design  and  the  data  produced 
by  each  observed  failure  is  comprised  of  the  4-tuple 

[j.  tj,  Siit^.Cj]  (26) 

representing  the  mode,  time,  stress  level,  and  design  configuration  of  the  failure. 

In  the  following  section,  four  approaches  cited  in  the  literature  for  analyzing 
accelerated  failure  time  data  are  discussed.  The  limitations  of  each  approach  are 
highlighted  with  regard  to  estimating  the  parameters  of  a  product’s  reliability 
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function  when  the  observed  data  is  of  the  form  shown  in  (26).  The  first  and  second 
subsections  describe  techniques  commonly  used  for  estimating  the  parameters  of  a 
product’s  reliability  function  where  a  specific  failure  mode  of  interest  has  been 
identified  prior  to  testing.  The  third  subsection  describes  an  approach  to  extend 
traditional  reliability  growth  models  to  incorporate  data  obtained  at  stress  levels 
that  are  not  representative  of  the  end  use  environment.  The  fourth  subsection 
discusses  a  commercial  database  used  to  predict  system  life  by  correlating  the  HALT 
results  of  multiple  products  with  their  subsequent  field  reliability.  Finally,  the  fifth 
subsection  suggests  a  direction  for  utilizing  data  of  the  form  shown  in  (26)  to 
improve  system  reliability  estimation  in  follow-on  testing. 

Modeling  Approaches 

Quantitative  Accelerated  Life  Testing  Model 

Quantitative  accelerated  life  testing  (QALT)  [45]  is  widely  used  to  expedite 
failures  and  quickly  acquire  information  regarding  a  product’s  time  to  failure 
distribution.  Test  planners  choose  the  type  and  levels  of  an  accelerating  variable 
along  with  the  fraction  of  samples  to  be  tested  at  each  level  to  minimize  the 
prediction  variance  of  product  life  with  respect  to  the  use-level  stress  (V-optimality) 
[60].  The  focus  of  QALT  is  typically  on  assessing  the  probability  of  failure  due  to 
specific  wearout  modes  since  a  majority  of  infant  mortality  failures  are  removed 
prior  to  QALT  through  product  ruggedization.  For  these  wearout  modes,  knowledge 
of  the  relationship  between  the  failure  mechanism  and  the  applied  stress  is  required 
and  may  be  available  through  an  understanding  of  the  chemical  or  physical 
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dynamics  causing  failure  or  through  previous  experience  with  similar  systems. 

Using  this  prior  information  enables  the  selection  of  an  appropriate  life-stress 
model  to  quantify  the  manner  in  which  the  time  to  failure  distribution  changes  with 
increased  stress. 

Analyzing  the  data  of  the  form  shown  in  (26)  using  QALT  models  requires 
coupling  an  appropriate  underlying  life  distribution  with  a  life-stress  relationship 
[45]  based  on  the  physical/chemical  dynamics  causing  failure  [54].  Any  distribution 
with  a  [0,  oo)  support  region  can  serve  as  the  underlying  life  distribution,  however 
selection  should  be  based  on  the  distribution’s  ability  to  accurately  describe  the 
data  as  shown  through  goodness  of  fit  testing  [61].  For  the  example  HALT  described 
above,  assume  that  the  underlying  life  distribution  relative  to  the  use-level  stress 
and  failure  mode  j  is  WEIB(a.j, )  with  CDF 


F(t;  ccj.pj)  =  P(T  <  t\a.j,Pj)  =  1  -  e  W  (27) 

where  the  value  of  the  shape  parameter  /?;  depends  only  on  the  failure  mode  j  [62], 
The  true  relationship  between  an  accelerating  variable  and  the  failure  mechanism 
can  often  be  quite  complicated,  thus  the  selection  of  a  life-stress  relationship  to 
model  a  failure  process  assumes  a  particular  failure  mode  of  interest  has  been 
identified  prior  to  testing.  Here,  we  assume  an  Inverse  Power-Law  stress-life 
relationship  [63,  64]  with  parameters,  K,n>  0  [55]  reflected  in  the  scale  parameter 
ocj.  Thus,  under  any  constant  stress  st  the  fraction  of  units  failing  due  to  mode  j 
prior  to  time  t  is 
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(28) 


Fj(t;  =  1  -  exp  {-(Kjs"Jt)Pl}. 


However,  in  a  step-stress  scenario  (28)  does  not  incorporate  the  exposure  across 
multiple  stress  levels  prior  to  failure  and  more  complex  techniques  [19,  65,  54]  are 
required.  Nelson  [18]  introduced  maximum  likelihood  estimation  to  step-stress 
QALT,  suggesting  a  Weibull-Inverse  Power  Law  likelihood  function 

x  fe-(*An;'((c-^))  7 

where  /  is  an  indicator  function  taking  on  a  value  of  one  if  the  observed  value  is  a 
failure  and  zero  it  is  a  censored  observation  and  incorporates  the  cumulative 
stress  exposure  prior  to  Sj.  Relating  failure  times  observed  at  elevated  stresses  to 
an  equivalent  use-stress  exposure  requires  at  least  three  observations  of  mode  j 
before  corrective  action  to  estimate  /?;,  Kj  and  n; .  But  this  requirement  violates  the 
iterative  test-fix-test  nature  of  HALT  which  only  allows  for  one  observation  for  each 
mode  prior  to  corrective  action.  Thus,  without  altering  the  HALT  philosophy,  the 
data  shown  in  (26)  is  insufficient  for  estimating  more  than  a  single  parameter  using 
any  life  stress  model  or  underlying  life  distribution. 


L  = 
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Physics  of  Failure  Analysis 


Physics  of  Failure  (PoF),  or  reliability  physics,  refers  to  a  scientific  approach 
wherein  modeling  and  simulation  are  used  to  understand  a  product’s  reaction  to 
external  stressors  and  address  its  fitness  for  use  with  respect  to  the  expected  use 
conditions  [66,  67].  The  intent  is  to  design-in  reliability  by  investigating  the 
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relationship  between  specific  failure  modes  and  accelerating  variables  and 
proactively  eliminating  the  root  causes  of  failure.  The  PoF  approach  provides  the 
strongest  characterization  of  system  reliability  [68],  incorporating  a  diverse  set  of 
engineering  disciplines  along  with  physics,  chemistry,  metallurgy,  mathematical 
statistics,  and  probability.  PoF  has  gained  wide  acceptance  among  military  and 
commercial  sectors,  showing  significant  time  and  cost  [69]  savings  while  allowing 
more  information  to  be  obtained  during  formal  testing. 

Pecht  et  al.  [70]  advocated  for  the  use  of  PoF  in  reliability  assessment  in  lieu  of 
the  popular  parts  count  technique  of  Mil  HDBK-217  [71].  Cushing  et  al.  [72] 
identified  several  limitations  of  the  parts  count  technique  that  could  be  addressed 
with  PoF  and  presented  procedures  for  implementing  the  approach  (Table  3). 
Mendel  [73]  introduced  probabilistic  physics  of  failure  (PPoF)  as  a  technique  to 
derive  the  statistical  lifetime  distribution  and  presented  a  case  for  applying  PPoF  in 
a  design  for  reliability  process.  Modarres  et  al.  [68]  further  emphasized  that  failure 
prediction  is  a  probabilistic  problem  due  to  uncertainties  associated  with  model 
parameters  and  failure-inducing  agents  that  can  result  from  changes  in 
environmental,  operating,  and  use  conditions. 

Alternatively,  Snook  et  al.  [74]  identifies  several  limitations  of  PoF,  most  notably 
that  for  immature  systems  or  those  with  multiple  usage  environments  PoF  analysis 
can  be  overly  complex  and  burdensome.  As  case  in  point,  Qi  [75]  presents  a  case 
study  evaluating  the  fatigue  life  of  printed  circuit  boards  under  thermal  cycling. 
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Table  3  -  Physics  of  failure  process 


Step 


Process 


5 

6 

7 


Define  realistic  requirements 

Define  the  expected  mechanical,  thermal,  electrical  and  chemical  usage 
environment  experienced  during  manufacture,  test,  operation,  storage,  and 
repair 

Identify  potential  failure  mechanisms  and  the  associated  degradation  processes 
(chemical,  electrical,  physical,  mechanical,  or  thermal] 

Define  appropriate  failure  models  and  input  parameters  relative  to  the  material 
characteristics,  damage  properties,  relevant  geometry,  and  operating 
environment 

Compute  the  variability  for  each  design  parameter  and  the  effective  reliability 
function. 

Design  to  the  usage  environment  incorporating  design  stress  spectra,  part  test 
spectra,  and  full-scale  test  spectra  based  on  the  anticipated  life-cycle 

Accept  the  design 


Thirty  thousand  ball  gate  array  (BGA)  interconnections,  arranged  in  154  packages, 
were  investigated  to  identify  which  were  likely  to  cause  the  circuit  boards  to  fall 
short  of  their  15-year  lifetime  requirement.  To  assess  the  likelihood  of  failure, 
three-dimensional  finite  element  models,  representing  each  distinct  BGA  package, 
were  exposed  to  fifteen  years  of  simulated  field  usage.  The  number  of  cycles  to 
failure  predicted  by  the  Engelmaier-Wild  solder  creep-fatigue  model  [76]  indicated 
that  multiple  BGA  packages  could  not  meet  the  fifteen  year  requirement. 

Qi  does  not  specify  what  corrective  action  was  implemented  subsequent  to  the 
investigation,  but  the  case  study  illustrates  that  the  level  of  effort  required  for  PoF 
necessitates  a  mature  design  and  a  narrowly  defined  failure  criteria.  Physics  of 
failure  simulations  are  tailored  to  replicate  specific  product  responses  under  a 
specific  set  of  usage  conditions  and  must  be  validated  against  demonstrated 
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performance.  Had  the  circuit  boards  contained  undiscovered  failure  modes,  the 
value  of  the  simulation  data  could  be  negated  if  corrective  action  alters  the  product’s 
response  to  the  applied  stress.  Identifying  these  latent  modes  through  simulation 
involves  replicating  product  responses  under  multiple  stress  environments  where 
validation  may  be  impossible  with  current  physical  models.  Thus  infant  mortality 
failures  must  be  removed  to  the  fullest  extent  possible  through  a  product 
ruggedization  test  such  as  HALT,  prior  to  starting  a  PoF  analysis. 

Accelerated  Reliability  Growth 

Reliability  Growth  is  the  positive  improvement  in  a  product’s  reliability  over 
time  due  to  design  changes  or  manufacturing  process  improvements  [9].  According 
to  Ebeling  [37],  the  earliest  developed  and  most  frequently  used  reliability  growth 
model  is  that  of  Duane  [38]  who  postulated  from  his  empirical  observations  that  for 
products  under  development,  plots  of  cumulative  failure  rate  versus  cumulative 
operating  hours  on  a  log-log  scale  were  approximately  linear  with  negative  slope. 
Crow  [39]  derived  the  stochastic  basis  for  the  Duane  model,  referred  to  in  the 
literature  as  a  power  law  process  (PLP)  or  informally  as  a  Weibull  process  -  a  type 
of  nonhomogeneous  Poisson  process  (NHPP)  with  a  Weibull  hazard  intensity 
function  h(t)  =  where  A  >  0  and  /?  >  0  are  interpreted  as  the  initial 

reliability  of  the  system  and  the  growth  rate  parameter,  respectively.  The  rigorous 
statistical  development  of  the  PLP  allows  estimation  and  prediction  in  both 
maximum  likelihood  and  Bayesian  contexts  and  has  made  the  approach  popular  in 
the  literature  [6,  37]. 
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Reliability  growth  testing  has  historically  existed  within  the  purview  of  large, 
complex  development  efforts,  such  as  military  weapon  systems  [77].  In  these 
programs  the  stress  profiles  used  in  testing  are  assumed  to  reflect  the  field-use 
environment.  Under  these  lower  stress  levels  testing  can  extend  several  months 
and  reliability  growth  planning  curves  (Figure  10)  are  used  to  ensure  adequate  test¬ 
time  is  made  available  to  discover  failures  and  subsequently  attain  the  reliability 
requirement  with  some  level  of  confidence.  Testing  exposes  several  distinct  failure 
modes,  thus  traditional  reliability  growth  models  track  only  the  overall  number  of 
failure  events  occurring  in  a  test  phase.  Reliability  improvement  is  therefore  shown 
as  a  decrease  in  the  rate  of  failure  events  across  test  phases  (assessed  growth).  An 
indication  that  the  reliability  growth  is  lagging  (assessments  falling  below  the 
planning  curve)  may  signal  the  need  to  reallocate  test  resources  or  reevaluate  the 
testing  strategy. 


Reliability 

(MTBF) 


Test  Phase  1  Test  Phase  2  Test  Phase  3 

Figure  10  -  Idealized  reliability  growth  planning  curves  developed  across  three  test 
phases  showing  the  difference  in  planned  and  assessed  reliability  (Source:  MIL- 

HDBK-189C,  p.  11) 
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In  commercial  environments  several  products  types  may  be  under  development 
simultaneously,  thus  little  time  may  exist  to  ruggedize  and  quantify  the  reliability  of 
a  particular  product.  In  these  situations,  accelerated  testing  methods  such  as  HALT 
may  be  utilized  as  part  of  the  reliability  improvement  process  prior  to  market 
introduction.  Like  traditional  growth  testing,  several  distinct  failure  modes  are 
observed,  but  many  of  these  observations  will  only  be  realized  at  elevated  stress 
levels.  The  Army  Materiel  Sustainment  and  Analysis  Activity  (AMSAA)  [9]  notes  that 
traditional  reliability  growth  models  produce  inaccurate  results  when  the  stress 
profile  is  not  operationally  relevant.  Strunz  and  Herrmann  [78]  address  planning, 
tracking  and  projecting  reliability  growth  of  liquid  rocket  engines  where  the  test 
profile  does  not  reflect  the  operational  use  environment.  Concluding  that  the 
traditional  modeling  approaches  are  insufficient,  the  authors  derive  a  Bayesian 
estimation  method  that  accounts  for  the  characteristics  of  the  test  profile  and 
aggregates  component,  subsystem  and  system  data. 

Feinberg  [77]  introduced  accelerated  reliability  growth  testing  (ARGT)  by 
modifying  the  AMSAA  [9]  reliability  growth  planning  model  equations  to  include  a 
system-level  acceleration  factor.  The  approach  is  based  on  the  assumption  of  a 
linear  correspondence  between  the  reliability  growth  attained  under  an  accelerated 
stress  and  that  which  would  have  occurred  using  only  the  use-level  stress.  Feinberg 
does  not  qualify  this  assumption,  implying  he  expects  it  to  hold  regardless  of  the 
intensity  of  the  accelerated  stress  and  the  level  of  extrapolation  considered. 
Intuitively,  as  stress  increases  the  time  to  failure  for  a  given  mode  decreases. 
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However,  as  noted  by  [45]  the  time  compression  is  not  equivalent  for  each  failure 
mode  in  the  system.  While  Feinberg  indicates  that  acceleration  factor  models  are 
mode  specific,  he  assumes  a  system-level  acceleration  factor  can  be  derived,  but 
leaves  this  task  to  the  reader. 

With  respect  to  estimating  product  reliability  from  the  data  in  (2),  reliability 
growth  models  improve  upon  QALT  and  PoF  in  that  design  changes  can  be 
accommodated.  However,  traditional  models  assume  growth  a  priori  via  a  fix 
effectiveness  value  applied  across  each  corrective  action.  In  an  accelerated  test,  this 
assumption  is  invalid  as  failures  may  be  discovered  far  from  the  use  stress.  Clearly, 
the  reliability  growth  attained  through  corrective  action  depends  on  the  likelihood 
of  the  failure  actually  occurring  in  operation  [54,  79, 17].  Nelson  [54]  discusses  this 
phenomenon,  describing  a  costly  effort  to  remove  a  failure  mode  which  never  would 
have  occurred  in  actual  use,  resulting  in  wasted  resources  and  no  actual  reliability 
improvement.  Moreover,  relating  growth  to  a  decrease  in  failure  events  is  also 
invalid  as  the  rate  of  failures  may  rise  with  increased  stress  even  though  reliability 
is  improving. 

McLean's  AFR  Estimator  Model 

After  collecting  HALT  data  from  more  than  fifty  products  across  twenty 
industries  McLean  [80]  developed  the  Actual  field  Failure  Rate  (AFR)  Estimator  to 
correlate  a  product’s  HALT  results  with  its  subsequent  field  failure  history.  When 
furnished  with  the  appropriate  HALT  data  and  product  information  (Table  4)  the 
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model  provides  the  field  failure  rate  for  the  improved  design  with  confidence  limits 
based  onj2  estimates  derived  from  the  SEMI  E10  standard  [81]. 


Table  4  -  AFR  Estimator  Inputs 


A 

B 

C 

D 

E 

F 

G 


Initial  MTBF  estimate  (early  testing,  reliability  prediction  standards) 
HALT  Hot  operating  limit  as  measured  on  the  product 
HALT  Cold  operating  limit  as  measured  on  the  product 
HALT  Vibration  operating  limit  as  measured  on  the  product 
Published  thermal  operating  specifications  (in  degrees  Celsius) 
Sample  size  used  in  the  final  HALT 
Field  duty  cycle 


The  AFR  Estimator  utilizes  linear,  quadratic  and  exponential  time  acceleration 
techniques  to  estimate  a  product’s  relative  life, 

MTBFredesign 

R  = - —  ('30') 

MTBForiginal  ^ 

representing  the  proportional  improvement  in  the  mean  time  between  failures  of 

the  redesigned  unit  over  that  of  the  original  design.  The  model  also  provides 

recommended  minimum  HALT  stress  levels  to  ensure  sufficient  test  coverage  in 

each  environment  and  to  assure  the  product  will  exceed  customer  expectations. 

Notwithstanding  these  advantages,  the  failure  rate  predictions  can  only  be  made  on 

the  basis  of  the  temperature  and  vibration  tests.  Further,  the  model  can't  estimate 

wear-out  mechanisms  which  must  be  addressed  using  QALT  techniques.  Although 

McLean  has  presented  results  [51]  of  the  AFR  Estimator,  the  model  remains 

unpublished  and  is  proprietary  to  Ops  A  La  Carte®  LLC.  Thus,  the  mechanics  of  how 
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the  linear,  quadratic  and  exponential  acceleration  techniques  are  utilized  is  unclear. 
Interestingly,  Brand  and  McLean  own  three  patents  [48,  49,  50]  outlining  models  for 
reliability  estimation  after  a  HALT  using  linear,  quadratic  and  exponential 
acceleration  techniques.  However,  correspondence  with  the  author  indicates  that 
the  AFR  Estimator  is  distinct  from  the  models  described  in  these  patents. 

Future  Research  Avenues 

Current  accelerated  reliability  techniques  investigate  the  response  of  a  known 
failure  mode  in  the  presence  of  one  or  more  accelerating  variables  assuming  the 
system  is  mature.  Estimating  reliability  from  data  obtained  from  qualitative 
accelerated  life  tests  [11,  42, 13]  depends  on  identifying  a  relationship  to  interpret 
the  observations  made  on  early  configurations  to  those  of  the  final  design.  This 
relationship  should  naturally  be  a  function  of  the  original  degradation  process,  but 
must  also  incorporate  changes  resulting  from  corrective  actions.  However,  for  even 
a  single  failure  mode-stressor  combination,  this  relationship  can  be  extremely 
complex  as  each  corrective  action  can  result  in  a  new  degradation  process  with 
unfamiliar  failure  states.  Although  reliability  physics  models  continue  to  improve, 
the  push  for  testing  chambers  that  can  produce  harsher  environments,  find  failures 
more  quickly,  and  test  more  complex  products  outpaces  the  capacity  to  derive  tools 
that  can  accurately  describe  product  behavior  in  these  environments. 

Consequently,  HALT  data  is  difficult  to  obtain  as  few  entities  have  the  equipment 
necessary  to  conduct  the  test  and  the  failure  information  is  often  considered 
competition  sensitive  to  both  the  item’s  manufacturer  and  the  testing  organization. 
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In  the  absence  of  sufficient  physical  models  and  available  data,  progress  in  this 
area  will  likely  be  confined  to  scanning  the  multi-dimensional  stress  space  bound  by 
the  HALT  destruct  limits  or  developing  novel  competing  risk  models  based  on 
reasonable,  but  unverifiable,  assumptions  on  the  form  of  the  data  and  model 
propriety  at  the  highest  stresses.  Escobar  and  Meeker  [60]  investigated  optimal 
accelerated  test  designs  when  two  or  more  independent  accelerating  variables  are 
present.  Gao  et  al.  [82]  compared  the  Escobar-Meeker  approach  to  orthogonal  and 
uniform  designs  in  computer  experiments  under  V-optimality.  Their  results  showed 
that  the  Escobar-Meeker  approach  performed  best  in  estimating  the  parameters  and 
the  pth  quantile  of  the  life  distribution  under  normal  stress. 

Assuming  independent  experimental  factors  may  be  reasonable  in  the 
neighborhood  of  the  design  stress.  But  at  limit  stress  failure  depends  on  the  total 
energy  put  into  the  system.  In  this  circumstance  the  design  space  is  non- 
quadrangular  with  an  outer  surface  representing  the  zero-sum  relationship 
between  the  factors.  Although  this  relationship  is  unknown  for  many  combinations 
of  experimental  factors,  comparing  the  Escobar-Meeker,  orthogonal  and  uniform 
designs  under  arbitrary  design  spaces  and  dependence  relationships  would  be  a 
valuable  contribution.  Additionally,  extending  methods  such  as  Calibrated 
Accelerated  Life  Testing  (CALT)  [47]  to  multiple  dimensions  may  prove  beneficial. 

Conclusion 

Maintaining  market  share  compels  manufacturers  to  provide  evolutionary  and 
revolutionary  capabilities  under  aggressive  development  schedules  while  ensuring 


56 


a  high  reliability  in  complex  stress  environments.  Successfully  attaining  reliability 
requirements  for  these  products  requires  setting  and  achieving  reliability  targets 
throughout  the  various  stages  of  the  development  process.  However,  schedule  and 
cost  constraints  often  preclude  sufficient  testing  on  early  prototypes  to  generate 
meaningful  reliability  estimates.  Under  these  constraints,  incorporating  data 
obtained  throughout  a  development  effort  can  result  in  improved  reliability 
estimation  prior  to  fielding.  Qualitative  accelerated  tests  are  often  used  to  quickly 
improve  system  reliability  by  identifying  and  removing  initial  design  flaws. 
Generally,  no  attempt  is  made  to  produce  a  reliability  measure  from  the  limited  data 
obtained  in  these  qualitative  tests  as  relevant  approaches  require  more  data.  The 
delineation  between  qualitative  accelerated  life  tests  and  QALT  exists  due  to  the 
philosophy  of  removing  the  root  cause  of  each  failure  through  corrective  action. 
However,  the  output  of  qualitative  testing  can  be  used  to  improve  reliability 
distribution  parameter  estimation  and  explore  the  outer  surface  of  the  enlarged 
design  space  to  guide  follow-on  tests. 
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IV.  Maximum  Likelihood  Estimation  for  the  Poly-Weibull  Distribution 
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ABSTRACT 

The  Weibull  distribution  has  long  been  a  popular  choice  for  modeling  lifetime  data 
of  various  mechanical  and  biological  phenomena  when  the  associated  hazard  rate 
function  is  constant  or  monotone  increasing  or  decreasing.  However,  non-monotone 
hazard  functions  are  common  in  reliability  and  survivability  contexts  where  a  system 
may  undergo  an  initial  “burn-in"  prior  to  periods  of  useful  life  and  eventual  wearout. 
In  these  scenarios,  the  standard  Weibull  can  only  model  a  portion  of  the  "bathtub” 
curve  but  is  incapable  of  adequately  modeling  the  entire  failure  process.  Several 
modifications  to  the  standard  two-parameter  Weibull  distribution  have  therefore 
been  introduced  in  the  literature  to  effectively  model  and  analyze  lifetime  data  where 
the  hazard  rate  function  is  bathtub  shaped.  The  performance  of  each  modified 
distribution  is  typically  assessed  by  its  ability  to  fit  reference  data  sets  that  are  known 
to  have  a  bathtub  shaped  hazard  rate  function.  The  current  paper  compares  the 
performance  of  two  recent  contributions  in  this  area  to  that  of  the  poly-Weibull 
distribution  with  respect  to  several  goodness  of  fit  measures.  In  addition,  numerical 
and  analytical  procedures  are  developed  for  obtaining  the  maximum  likelihood 
parameter  estimates  and  standard  errors  for  the  generalized  poly-Weibull 
distribution  with  arbitrary  number  of  terms.  Our  results  show  that  both  the  bi- 
Weibull  and  tri-Weibull  distributions  fit  the  reference  data  sets  better  that  either  of 
the  current  best-fit  models. 


Introduction 

The  Weibull  distribution  has  long  been  a  popular  choice  for  modeling  lifetime 
data  of  various  mechanical  and  biological  phenomena  when  the  associated  hazard 
rate  function  is  either  constant  or  monotone  increasing  or  decreasing.  However, 
non-monotone  hazard  functions  are  common  in  reliability  and  survivability  contexts 
where  a  system  may  undergo  an  initial  "burn-in"  prior  to  periods  of  useful  life  and 
eventual  wearout.  In  these  scenarios,  the  standard  Weibull  distribution  can  model  a 
portion  of  the  "bathtub”  hazard  curve  but  is  incapable  of  describing  the  entire 
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failure  process.  Several  modifications  to  the  standard  two-parameter  Weibull 
distribution  have  therefore  been  introduced  in  the  literature  to  effectively  model 
lifetime  data  where  the  hazard  rate  function  is  bathtub  shaped.  The  performance  of 
each  modification  is  assessed  by  comparing  its  goodness  of  fit  to  data  sets  published 
by  Aarset  [83]  and  Meeker  [84],  known  to  have  bathtub  shaped  hazard  functions. 

Bathtub-shaped  hazard  functions  arise  from  the  existence  of  multiple  competing 
failure  modes  which  dominate  at  different  epochs  in  the  life  of  a  system.  The  poly- 
Weibull  [85]  distribution  arises  naturally  in  scenarios  of  competing  risks  as  it 
describes  the  minimum  of  several  independent  random  variables  where  each 
follows  a  distinct  Weibull  law.  The  current  paper  compares  the  goodness  of  fit  of 
the  poly-Weibull  to  two  recently  proposed  Weibull  modifications,  the  new  modified 
Weibull  distribution  [86]  and  the  exponentiated  modified  Weibull  extension 
distribution  [87].  Our  results  show  that  the  poly-Weibull  fits  both  data  sets  better 
than  either  of  these  modifications. 

Almalki  and  Yuan  [86]  introduced  the  new  modified  Weibull  (NMW)  distribution 

FNMW(t)a’8,P,Y’h)  =  1  -  e~ate~PtYeM  a,  9,  /?,  y,  A  >  0,  t>0  (31) 

by  considering  a  two-component  serial  arrangement  in  which  one  component 
follows  a  standard  two-parameter  Weibull  model  and  the  other  follows  a  modified 
Weibull  (MW)  distribution  [88].  The  authors  show  that  the  goodness  of  fit  of  the 
NMW  is  superior  to  that  of  other  four  and  five  parameter  Weibull  modifications 
including  the  additive  Weibull  distribution  [89],  the  modified  Weibull  distribution 


59 


[90]  and  the  beta  modified  Weibull  distribution  [91].  Likewise,  Sarhan  and  Apaloo 
[87]  showed  that  the  exponentiated  modified  Weibull  extension  (EWME)  model 

FEWME(.t'-A-’a’P’Y)=  1  -  eAa(1_e<  /  >  )  A,a,p,Y>  0,  t>  0  (32) 

fit  the  reference  data  better  than  the  exponentiated  Weibull  [92],  the  exponentiated 
Gompertz  [93]  and  the  modified  Weibull  extension  distributions  [94]. 


The  Poly-Weibull  Distribution 

The  CDF  of  the  poly-Weibull  distribution  is  expressed  as: 


FPW(t\  «,/?)  =  1  —  |  exp  -  ^  i  <Xj, >  0,  t  >  0 


where  ctj,  /?;  represent  the  scale  and  shape  parameters  associated  with  the  Weibull 
model  describing  risk  j  —  1, ... ,].  Accordingly,  the  poly-Weibull  density  function  is 

fPW(t\a,p)  =  R(t)h(t)  =  jexp  j  *  ^ ^  p.  ■  (34) 

When  /  =  2,  equations  (33)  and  (34)  are  the  CDF  and  pdf  of  the  bi-Weibull 
distribution,  and  when  J  —  3  the  model  is  naturally  known  as  the  tri-Weibull 
distribution.  Considering  the  poly-Weibull  was  introduced  over  twenty  years  ago 
[85],  the  body  of  literature  studying  its  properties  is  quite  limited.  Within  this  body, 
the  predominant  focus  has  been  on  Bayesian  estimation  of  the  bi-Weibull  model 
parameters.  Berger  and  Sun  [85]  studied  the  bi-Weibull  distribution  using  a  Gibbs 
sampling  algorithm  with  informative  and  non-informative  priors.  The  authors 
computed  the  posterior  density  and  predictive  reliability  when  the  shape 


parameters  are  either  known  or  unknown.  Their  simulation  results  compared 
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favorably  with  a  closed-form  exact  solution,  but  were  shown  to  be  computationally 
expensive  for  small  sample  sizes.  Davison  and  Louzada-Neto  [95]  also  explored  the 
bi-Weibull  model,  however,  they  assert  that  Markov  chain  simulation  is  not 
necessary  for  obtaining  posterior  probabilities.  Using  real  and  generated  data,  the 
authors  illustrate,  that  either  Laplace’s  method  [96]  or  the  Bayesian  bootstrap  [97] 
are  sufficient  and  can  reduce  the  computational  burden. 

Figures  11  and  12  illustrate  several  bathtub  shapes  modeled  with  the  bi-Weibull 
and  tri-Weibull  distributions  with  parameter  vector  6  —  [/31, ... ,  /?y,  av  ... ,  a.j]. 


Figure  11  -  Example  bi-Weibull  density  (top)  and  hazard  (bottom)  functions 
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Comparing  the  center  sections  of  each  figure,  the  reader  will  note  the  ability  of  the 
tri-Weibull  to  model  complex  failure  processes  for  which  the  "useful  life"  portion  of 
the  hazard  function  is  not  constant  and  thus  does  not  follow  an  exponential  model. 


Figure  12  -  Example  tri-Weibull  density  (top)  and  hazard  (bottom)  functions 


An  advantage  of  the  poly-Weibull  model  is  that  the  parameter  values  have  the  same 
relationship  to  the  mean  and  variance  of  the  failure  process  as  does  the  standard 
two-parameter  Weibull  model.  Thus  <  1  implies  a  failure  process  with  a  large 
coefficient  of  variation  and  decreasing  hazard  rate  indicative  of  infant  mortality. 
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Conversely,  /?;  >  1  portends  a  wear-out  failure  mechanism  with  an  increasing 
hazard  rate  function  reflective  of  a  lower  coefficient  of  variation.  The  poly-Weibull 
distribution  is  therefore  capable  of  generating  bathtub  shaped  hazard  functions  by 
modeling  multiple  failure  processes  simultaneously  as  Figures  11  and  12  illustrate. 


Parameter  Estimation 

Consider  a  random  sample  of  observations  t1(  ...,tn  from  a  poly-Weibull(a,/?) 
distribution  with  unknown  parameter  vector  6  —  (a1; ... ,  a]t  /31, ... ,  /?y),  and  indicator 
variable  8t  where  8t  —  1  if  k  is  a  failure  time  and  8t  —  0  if  tj  is  a  censoring  time.  The 
log-likelihood  function  £(0)  based  on  data  (tv  8J, ... ,  (tn,  8n )  is 


i=i 


J  't-\^ 

Ll 


1(6)  =  £  5,  log^ 


]= 1 


;= i 


(35) 


Davison  and  Louzada-Neto  [95]  initially  presented  the  system  of  2J  nonlinear 
equations,  which  when  solved  would  provide  the  maximum  likelihood  estimates  for 
ccj  and  j3j,  j  —  1, ... ,].  However,  these  equations  contained  an  error  -  omitting  /?; 
from  the  denominator  of  dL/daj.  The  corrected  equations  for  obtaining  the  poly- 
Weibull  MLE’s,  which  have  been  verified  both  analytically  and  numerically,  are 
therefore 
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where  hj{t{)  =  p.t^i  1ctj  ^  and  h(t{)  =  T/j=1  (fij t?J  1a.  ^  j.  Solving  this  system  of  non¬ 
linear  equations  cannot  be  accomplished  analytically  and  the  use  of  Newtonian  or 
quasi-Newtonian  numerical  optimization  techniques  can  be  tedious  as  finding  a 
solution  is  highly  sensitive  on  the  starting  values  for  the  parameters  in  each 
equation.  We  find  it  simpler  to  obtain  accurate  parameter  estimates  by  maximizing 
the  log-likelihood  function  directly  using  a  quasi-Newtonian  algorithm  [98]. 
However,  for  asymptotic  interval  estimation,  the  optimization  algorithm  can 
produce  inaccurate  Hessian  matrices  leading  to  negative  values  along  the  diagonal 
of  the  covariance  matrix.  Thus  for  finding  the  standard  errors  of  the  poly-Weibull 
model  parameters  the  components  of  the  observed  Fisher  information  matrix 
0(9*)  —  — VVT£(0)|0=e*  have  been  derived  analytically  and  are  presented  in 
appendix  A. 

Application 

In  this  section  the  Aarset  [83]  and  Meeker  [84]  data  sets  are  analyzed  and  the 
goodness  of  fits  for  the  bi-Weibull  and  tri-Weibull  distributions  are  compared  to 
those  of  the  NMW  [86]  and  the  EMWE  [87].  In  keeping  with  the  established 
precedent,  the  goodness  of  fit  comparison  in  this  paper  is  based  on  the  values  of  the 
Kolmogorov-Smirnov  (K-S)  test  statistic,  the  Akaike  information  criterion  (AIC)  and 
the  log-likelihood  function  for  each  model  computed  at  their  respective  MLE's.  In 
addition,  the  fit  of  each  distribution  is  assessed  graphically  as  the  reliability,  density 
and  hazard  functions  are  plotted  against  their  nonparametric  counterparts. 
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Aarset  data 


The  Aarset  [83]  data  set  (Table  5)  represents  the  lifetimes  of  50  devices  and 
contains  no  censored  observations. 

Table  5  -  The  Aarset  data  set 
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Table  6  shows  the  MLE’s  and  standard  errors  for  the  parameters  of  each  of  the  four 
models  considered  while  Table  7  displays  a  comparison  of  each  model’s  goodness  of 
fit  measures.  The  data  in  Table  7  indicates  that  the  null  hypothesis  of  the  two- 
sample  K-S  test  cannot  be  rejected  for  any  of  the  four  models  at  a  significance  level 
below  0.8.  However,  the  data  also  shows  that  the  tri-Weibull  and  bi-Weibull  fit  the 
data  better  than  either  the  NMW  or  the  EMWE  as  both  have  larger  likelihoods  as 
well  as  smaller  K-S  statistics  and  AIC  values.  Further,  the  superior  performances  of 
the  poly-Weibull  models  are  immediately  clear  upon  observing  Figure  13.  In  the  top 
plot,  the  reliability  function  of  each  model  is  plotted  against  the  Kaplan-Meier  [30] 
nonparametric  estimate  of  the  reliability  function  for  the  data. 
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Table  6  -  Maximum  likelihood  estimates  (standard  errors)  for  the  Aarset  data 


Model 

MLE  of  the  Parameters 

k 

k 

k 

a-L 

a2 

a3 

Tri-Weibull 

98.152 

0.524 

4.215 

85.091 

122.478 

92.299 

[32.762] 

[0.056] 

[0.937] 

[0.339] 

[52.740] 

[8.870] 

k 

k 

a-L 

a  2 

Bi-Weibull 

82.334 

0.702 

84.907 

61.663 

[22.602] 

[0.075] 

[0.328] 

[14.538] 

a 

X 

P 

Y 

EMWE 

49.050 

7.18xl0-5 

3.148 

0.145 

* 

* 

* 

* 

a 

X 

P 

Y 

e 

NMW 

0.071 

0.197 

7.015xl0-8 

0.016 

0.595 

[0.031] 

[0.184] 

[1.501x10-7] 

[3.602] 

[0.128] 

The  NMW  and  EWME  fit  the  central  portion  of  the  data  better  that  the  bi-Weibull 


while  both  the  bi-Weibull  and  tri-Weibull  more  closely  fit  the  upper  and  lower  tails 


of  the  data  where  the  majority  of  observations  are  concentrated.  However,  the  tri- 


Weibull  is  clearly  the  best  fit  throughout  the  entire  range  of  the  observations. 


Table  7  -  Performance  measures  for  the  Aarset  data  set 


Model 

Params 

Log-Lik 

K-S 

p-value 

AIC 

Tri-Weibull 

6 

-202.51 

0.063 

0.998 

417.01 

Bi-Weibull 

4 

-206.09 

0.100 

0.925 

420.20 

EMWE 

4 

-213.86 

0.101 

0.646 

435.72 

NMW 

5 

-212.90 

0.088 

0.803 

435.80 

In  the  middle  plot,  each  model’s  density  function  is  plotted  against  a  histogram  of 
the  data.  Note  that  the  poly-Weibull  models  indicate  that  the  probability  of  failure 
after  the  final  observation  is  near  zero,  as  would  be  expected  for  a  system  with  a 
true  bathtub-shaped  hazard,  while  the  NMW  and  EMWE  do  not  reflect  this.  This  is 
also  clearly  evident  in  the  bottom  plot  where  the  hazard  functions  are  plotted 
against  the  empirical  hazard  plot  [99]. 
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Figure  13  -  Triple  plot  illustrating  the  fit  of  the  bi-Weibull,  tri-Weibull  NMW 
and  EMWE  models  to  the  Aarset  data 
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Meeker  data  -  Uncensored 


The  Meeker  [84]  data  set  (Table  8)  represents  the  observed  lifetimes  of  30 
devices  and  includes  eight  censored  observations.  The  data  set  only  has  a  bathtub 
shaped  hazard  function  if  the  censored  observations  are  treated  as  failures.  Thus, 
the  subsequent  analysis  follows  this  precedent  established  in  the  literature.  In  the 
following  subsection  the  Meeker  data  is  analyzed  again,  correctly  treating  the  final 
observations  as  censored. 


Table  8  -  The  Meeker  data  set 
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Maximum  likelihood  parameter  estimates  and  standard  errors  are  displayed  in 
Table  9.  Note  the  large  standard  errors  associated  with  the  bi-Weibull  and  tri- 
Weibull  shape  parameter  /?;=1.  As  a  result  of  treating  the  censored  observations  in 
the  Meeker  data  set  as  failures  we  have  enforced  a  deterministic  failure  at  300  hours 
such  that 

P(T  >  300|T  =  300)  =  0.  (37) 

Thus,  the  optimization  procedure  attempts  to  fit  bi-Weibull  and  tri-Weibull 
distributions  to  the  uncensored  Meeker  data  with  shape  parameter  /?x  ->  oo.  The 
large  standard  errors  therefore  reflect  the  difference  between  the  true  distribution 
governing  the  failure  process  and  the  distribution  that  can  be  numerically  obtained. 
Observing  the  density  curves  in  the  center  plot  of  Figure  14,  shows  that  the  poly- 
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Weibull  accurately  portrays  a  deterministic  failure  at  300  hours  while  the  NMW  and 


EWME  curves  incorrectly  show  that  failures  beyond  300  hours  are  possible. 


Table  9  -  Maximum  likelihood  estimates  (standard  errors)  for  the  uncensored 

Meeker  data  set 


Model 

MLE  of  the  Parameters 

k 

k 

k 

at 

a2 

«3 

Tri-Weibull 

124.000 

5.591 

0.738 

299.850 

356.022 

352.185 

(188.557) 

(1.907) 

(0.102) 

(1.825) 

(46.686) 

(136.733) 

Bi-Weibull 

k 

124.000 

k 

0.892 

a1 

299.780 

a2 

253.940 

(98.228) 

(0.122) 

(1.117) 

(62.099) 

EMWE 

a 

197.210 

X 

5.468  xlO-6 

0 

4.482 

Y 

0.129 

* 

* 

* 

* 

a 

X 

0 

Y 

6 

NMW 

0.024 

0.056 

5.991  xlO-8 

0.012 

0.629 

(0.019) 

(0.024) 

(8.164x10-8) 

(1.290) 

(0.158) 

Regardless,  the  goodness  of  fit  measures  in  Table  10  indicate  that  bi-  and  tri-Weibull 


remain  the  preferred  models  when  compared  to  the  NMW  and  EMWE  for  fitting  the 


data.  With  additional  computing  power  this  fit  can  be  improved. 


Table  10  -  Performance  measures  for  the  uncensored  Meeker  data  set 


Model 

Params 

Log-Lik 

K-S 

p-value 

AIC 

Bi-Weibull 

4 

-156.14 

0.105 

1.000 

320.30 

Tri-  Weibull 

6 

-154.65 

0.097 

1.000 

321.30 

EMWE 

4 

-166.35 

0.131 

0.632 

340.71 

NMW 

5 

-166.18 

0.148 

0.482 

344.40 

The  graphical  fit,  as  shown  in  Figure  14,  again  indicates  that  both  the  tri-Weibull 
and  the  bi-Weibull  outperform  the  NMW  and  EMWE.  However,  because  the  bathtub 
shape  is  less  marked  for  the  uncensored  Meeker  data,  the  difference  between  the  bi- 
and  tri-Weibull  is  insufficient  to  overcome  the  penalty  of  two  additional  parameters. 
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Figure  14  -  Triple  plot  illustrating  the  fit  of  the  bi-Weibull,  tri-Weibull  NMW 
and  EMWE  models  to  the  uncensored  Meeker  data 


70 


Meeker  data  -  Censored 


The  Meeker  [84]  data  set  is  analyzed  again,  this  time  treating  the  final  eight 
observations  as  right  censored.  As  the  respective  authors  do  not  present  results  for 
the  NMW  and  EMWE  in  this  scenario,  the  results  presented  below  were  obtained 
using  the  numerical  procedure  discussed  above  for  the  poly-Weibull  distribution. 
The  maximum  likelihood  parameter  estimates  are  shown  in  Table  11. 


Table  11  -  Maximum  likelihood  estimates  for  the  censored  Meeker  data  set 


Model 

MLE  of  the  Parameters 

Tri-Weibull 

k 

6.795 

Pz 

0.742 

@3 

22.125 

ax 

338.686 

a2 

346.727 

«3 

974.718 

Bi-Weibull 

k 

6.795 

$2 

0.742 

a1 

338.686 

a2 

346.727 

EMWE 

a 

93.100 

X 

2.2x10-3 

P 

0.633 

Y 

0.625 

NMW 

a 

0.0142 

X 

0.0346 

P 

0.6939 

Y 

0.3015 

6 

0.0117 

Here,  /?1;  /?2,  cLlt  a2  are  equivalent  for  the  bi-Weibull  and  tri-Weibull  distributions  as 


/?3  and  a3  provide  no  additional  information,  resulting  in  a  singular  Hessian  for  the 
tri-Weibull.  Thus,  the  goodness  of  fit  measures  displayed  in  Table  12  show  that  bi- 
Weibull  is  the  preferred  model. 

Table  12  -  Performance  measures  for  the  censored  Meeker  data  set 


Model 

Params 

Log-Lik 

K-S 

p-value 

AIC 

Bi-Weibull 

4 

-140.95 

0.1304 

0.9924 

289.9 

NMW 

5 

-141.16 

0.1304 

0.9924 

292.5 

Tri-Weibull 

6 

-140.95 

0.1304 

0.9924 

293.9 

EMWE 

4 

-151.43 

0.2609 

0.4218 

310.9 

Figure  15  shows  that  the  poly-Weibull  and  NMW  fit  the  data  equally  well,  while  the 
performance  of  the  EMWE  reflects  convergence  problems  associated  with  the 
model.  Note  that  in  each  plot  the  bi-  and  tri-Weibull  curves  are  coincident. 
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Figure  15  -  Triple  plot  illustrating  the  fit  of  the  bi-Weibull,  tri-Weibull  NMW 
and  EMWE  models  to  the  censored  Meeker  data 
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Conclusion 


The  poly-Weibull  distribution  was  presented  as  an  attractive  alternative  to  the 
standard  two-parameter  Weibull  for  fitting  data  sets  where  the  hazard  function  is 
monotone,  non-monotone  or  even  bathtub  shaped.  Numerical  and  analytical 
procedures  for  obtaining  the  poly-Weibull  maximum  likelihood  parameter 
estimates  and  asymptotic  standard  errors  are  presented.  It  was  shown  that  the 
equations  for  obtaining  the  poly-Weibull  maximum  likelihood  parameter  estimates 
initially  derived  and  presented  by  Davison  and  Louzada-Neto  contained  an  error. 
The  corrected  equations  were  therefore  derived  and  verified  along  with  the 
components  of  the  observed  Fisher  Information  matrix  for  the  generalized  poly- 
Weibull  distribution.  The  goodness  of  fit  for  two  forms  of  the  poly-Weibull 
distribution,  the  bi-Weibull  and  tri-Weibull,  was  then  assessed  for  two  reference 
data  sets  using  the  Akaike  information  criterion  and  Kolmogorov-Smirnov  test 
statistic.  Our  results  show  that  the  bi-Weibull  and  the  tri-Weibull  outperform  two 
other  modified  Weibull  distributions  with  respect  to  their  fit  of  data  for  which  the 
hazard  rate  function  is  bathtub  shaped.  Further,  the  tri-Weibull  dominates  the  other 
models  when  the  bathtub  shape  is  prominent.  However,  this  advantage  is  overcome 
by  the  penalty  incurred  in  the  AIC  from  additional  parameters  when  the  bathtub 
shape  is  less  apparent.  Thus  it  is  our  recommendation  that  the  tri-Weibull 
distribution  be  considered  when  one  needs  to  model  data  with  a  pronounced 
bathtub  shaped  hazard  function.  When  the  data  does  not  reflect  a  prominent 
bathtub  shape  the  bi-Weibull  distribution  is  an  excellent  choice. 
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V.  Characterizing  Reliability  Growth  in  Early  System  Design 


Jason  K.  Freels 
Joseph  J.  Pignatiello 
Richard  L.  Warr 
Raymond  R.  Hill 


ABSTRACT 

A  successful  reliability  growth  program  requires  setting  and  achieving  reliability 
targets  throughout  the  various  stages  of  a  system’s  development  process.  However, 
schedule  and  cost  constraints  often  preclude  sufficient  testing  on  early  prototypes  to 
generate  meaningful  reliability  estimates.  While  qualitative  accelerated  test  methods 
quickly  improve  system  reliability  by  identifying  and  removing  initial  design  flaws, 
no  attempt  is  made  at  quantifying  the  reliability.  In  the  current  paper  a  modified 
accelerated  life  test  is  proposed  whereby  the  projected  system  reliability  can  be 
estimated  after  implementing  corrective  action.  Assuming  the  system  contains  an 
unknown  number  of  independent  competing  failure  modes  whose  respective  time  to 
occurrence  is  governed  by  a  distinct  Weibull  law,  the  observed  failure  times  are 
modeled  with  the  poly-Weibull  distribution.  We  show  that  under  a  qualitative 
accelerated  life  test  scenario  the  poly-Weibull  failure  process  can  be  modeled  with  a 
Weibull  distribution  for  various  sample  sizes  and  system  types.  Thus,  the  proposed 
model  utilizes  the  Weibull  distribution  to  estimate  system  reliability  after  one  or 
more  failure  modes  have  been  discovered  and  removed.  To  our  knowledge  the 
proposed  model  is  the  first  attempt  to  incorporate  physical  acceleration  models  into 
reliability  growth  planning  and  is  intended  to  serve  as  a  prototype  upon  which 
refinements  in  reliability  growth  modeling  can  be  incorporated. 

Introduction 

Reliability  growth  testing  has  typically  existed  within  the  purview  of  large 
development  efforts,  such  as  military  weapon  systems  [77]  where  failures  may  not 
be  observed  until  the  system  is  inspected  at  the  end  of  a  testing  phase.  In  this 
context,  the  metric  of  interest  for  many  reliability  growth  models,  such  as  those 
based  on  the  power  law  process  (PLP)  [39],  is  the  cumulative  mean  time  between 
failures  within  a  given  phase  [9].  Accordingly,  these  models  express  reliability 
growth  as  a  change  in  the  rate  of  occurrence  of  failures  (ROCOF)  across  subsequent 
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test  phases.  But  reliability  growth  testing  can  take  several  months  or  longer  to 
reach  a  reliability  goal  with  a  specified  level  of  confidence  as  PLP-based  models 
assume  the  test  conditions  reflect  the  intended  usage  environment.  Consequently, 
accelerated  testing  methods  have  been  viewed  as  a  means  to  reduce  the  duration  of 
a  reliability  growth  test. 

Introducing  accelerated  stresses  in  a  reliability  growth  test  can  increase  the 
ROCOF  as  the  energy  put  into  the  system  is  sufficient  to  activate  failure  modes  that 
would  otherwise  remain  dormant  under  use-level  stress.  However,  no  model  exists 
to  describe  how  changes  in  the  level  of  an  accelerating  variable  affect  the  ROCOF. 
Moreover,  when  multiple  failure  modes  are  exposed  in  a  test,  separate  models  are 
required  to  characterize  the  life-stress  relationship  for  each  mode  [54]  as  models 
describing  the  behavior  of  entire  systems  under  accelerated  stress  also  do  not  exist. 

One  particular  testing  strategy,  qualitative  accelerated  reliability  testing  [12, 11, 
13],  exposes  early  product  designs  to  elevated  environmental  and  repetitive  use 
stresses  to  force  latent  defects  to  become  manifest  thereby  ensuring  that  with 
appropriate  mitigation  a  reliable  product  is  fielded  quickly.  However,  qualitative 
test  methods  are  designed  to  discover  and  remove  failure  modes  quickly,  without 
regard  for  reliability  estimation  and  many  tests  employ  a  test-fix-test  strategy 
where  only  a  single  observation  of  each  failure  mode  is  obtained  prior  to  its  removal 
through  corrective  action.  Test  planners  have  long  sought  a  way  to  incorporate  the 
results  of  qualitative  accelerated  testing  into  an  estimate  of  system  reliability  in  the 
early  design  and  development  phases.  While  case-studies  [15, 16, 14]  attest  to  the 
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effectiveness  of  qualitative  testing  for  producing  reliable  products,  translating  the 
limited  failure  data  into  a  meaningful  measure  of  reliability  growth  has  remained 
elusive.  Unlike  qualitative  testing  which  discovers  previously  unknown  failure 
modes  for  which  little  information  on  the  mode-specific  time  to  failure  distributions 
may  be  available,  quantitative  accelerated  life  tests  (QALT)  investigates  the 
probability  of  failure  due  to  a  few  known  failure  modes. 

In  a  quantitative  test,  prior  knowledge  of  the  relationship  between  the  failure 
mechanism  and  the  accelerating  variable  is  required  to  identify  an  appropriate  life- 
stress  model  used  to  describe  how  the  time  to  failure  distribution  changes  with 
increased  stress  [54,  45].  For  qualitative  tests,  prior  knowledge  is  unavailable  since 
the  failure  modes  may  be  unknown  prior  to  their  occurrence  and  the  stresses  used 
to  discover  the  failures  often  follow  complex  profiles  [17].  Further,  corrective 
action  alters  the  system  design  configuration  by  eliminating  the  flaws  causing 
failure.  Since  subsequent  configurations  represent  different  systems,  their 
respective  failure  data  cannot  be  combined  in  a  meaningful  way.  To  address  this 
issue  the  current  paper  proposes  a  modified  qualitative  test  and  an  associated 
model  whereby  observations  obtained  in  an  accelerated  test  can  be  used  to  estimate 
the  projected  system  reliability  after  corrective  action. 

Testing  Framework 

In  the  proposed  testing  scenario  a  sample  of  n  prototypes,  sharing  a  common 
initial  design  configuration  are  exposed  to  a  constant  elevated  stress  state  until 
failure  (i.e.,  there  is  no  censoring).  The  stress  state  may  be  comprised  of  one  or 
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more  stressors,  but  the  levels  of  each  stressor  must  remain  constant  throughout  the 
test  sequence  until  all  n  prototypes  have  failed.  Once  the  test  is  complete,  corrective 
action  is  implemented  on  all  of  the  prototypes  in  the  sample  set  to  remove  each 
distinct  failure  mode  discovered  thereby  restoring  the  prototypes  to  common 
design  configuration.  In  subsequent  tests,  the  levels  of  one  or  more  stressors  may 
be  changed  to  create  a  more  extreme  stress  state  capable  of  generating  additional 
failures.  This  process  continues  until  a  failure  occurs  that  is  deemed  too  cost 
prohibitive  to  remove  at  which  point  the  test  is  concluded  and  the  design 
configuration  is  frozen.  In  the  next  section  we  introduce  the  model,  first  discussing 
the  analysis  of  data  from  an  accelerated  test  in  which  each  system  contains  a  single 
failure  mode  and  no  corrective  action  is  implemented.  Next,  the  data  are  re¬ 
analyzed  assuming  that  corrective  action  is  implemented  after  the  test.  Finally,  we 
consider  the  general  case  where  multiple  failure  modes  exist  in  each  system  and 
corrective  action  is  implemented. 

Model 

Single  Failure  Mode  Case  -  No  Corrective  Action 

In  the  simplest  case,  only  one  failure  mode  exists  in  each  system.  For  this  case 
the  n  times-to-failure  can  be  modeled  by  the  distribution  function  F(t;  6 )  where  the 
parameter  vector  is  a  function  of  stress,  0(s).  In  the  analysis  of  accelerated  test  data 
it  is  often  assumed  that  the  relative  standard  deviation  (er  ■  /U-1)  of  the  failure 
observations  is  unaffected  by  changes  in  the  applied  stress  [100].  For  shape-scale 
distributions  such  as  the  Weibull(a,  /?),  gamma(a,  /?),  and  log-normal(eM,  a2)  this 
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implies  that  the  functional  relationship  between  a  product’s  lifetime  distribution 
and  the  applied  stress  is  reflected  in  the  scale  parameter  only,  while  the  shape 
parameter  is  independent  of  stress. 

When  the  failures  from  an  accelerated  life  test  are  due  to  a  single  failure  mode 
the  observations  can  be  rank  ordered  from  smallest  to  largest  and  plotted  against 
their  corresponding  quantiles  as  derived  from  a  nonparametric  estimate  of  the 
cumulative  distribution  function,  F(t).  This  nonparametric  estimate  can  be  found 
using  a  number  of  plotting  position  models  [101,  pp.  6-8],  such  as  the  median 
plotting  position 

(i  —  0.3  \ 

F(t‘)  =  100  ■  tr 04  )  (38) 

where  i  denotes  the  ith  ordered  failure  and  n  represents  the  total  number  of 
observations.  To  determine  if  the  data  follow  a  specific  candidate  distribution,  the 
plotting  coordinate  values  and  failure  times  are  substituted  into  the  "linearized" 
distribution  function.  If  the  data  are  derived  from  the  candidate  model  the 
transformed  observations  should  closely  trace  a  line.  Using  the  Weibull(a,  /?) 
distribution  as  an  example,  the  linearized  distribution  function  is  found  by  rewriting 
the  CDF  as 


logio 


=  /?  logic (ti)  -/?log10 (a). 


(39) 


If  the  data  are  consistent  with  the  Weibull  model,  the  resulting  plot  of 
—  In  (l  —  F(tj)^  versus  tt  will  be  nearly  linear  with  slope  /?  when  plotted  on  a  log- 
log  scale  (Figure  16).  Alternatively,  one  may  also  use  specially  designed  plotting 
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papers  that  have  been  developed  for  several  distributions  where  the  axes  have  been 
transformed  such  that  the  untransformed  failure  times  and  non-parametric 
probability  estimates  F(tj)  may  be  plotted  directly.  If  Weibull  plotting  paper  is  used 
a  and  /?  may  be  estimated  graphically,  observing  where  the  fit  line  crosses  the  line 
designating  the  0.632  quantile  and  the  slope  of  the  plot.  Figure  16  displays  the  raw 
data  and  linear  fits  from  an  accelerated  life  test  discussed  in  [54,  p.  115]  where  the 
results  are  hours  to  failure  obtained  at  three  levels  of  temperature.  Because  /?220  = 
/?24o  —  /?26o  —  P  was  assumed,  the  fit  lines  should  be  nearly  parallel.  Kececioglu 
[55]  notes  that  for  small  data  sets  (i.  e. ,  those  where  n  <  10  samples)  the  slope  of 
the  distribution  fit  lines  can  vary  widely  from  their  true  values.  Nelson  [54] 
suggests  that  unless  the  line  slopes  change  systematically  with  stress  or  the  slope 
for  one  stress  level  is  dramatically  different  than  the  others,  the  lines  may  be 
redrawn  parallel  to  fit  with  the  constant  /?  assumption.  From  Figure  16  it  is  clear 
that  the  slopes  of  the  observations  obtained  at  220°C  and  240°  C  are  nearly  parallel 
and  the  slope  of  the  260°  data  is  not  dramatically  different.  Thus  the  260 °C  fit  line 
is  redrawn  such  that  all  three  fit  lines  are  parallel  and  the  median  time  to  failure  is 
unchanged  as  shown  in  Figure  17.  From  these  adjusted  plots  it  is  determined  that 
/?  -  6  and  a220  —  2850,  a240  =  1700,  a260  =  1040. 
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Figure  16  -  Weibull  plot  of  failure  times  observed  at  three  levels  of 
temperature  in  an  example  accelerated  life  test  [18] 


Probability,  F(t) 

Figure  17  -  Weibull  plot  of  example  failure  data  [18]  redrawn  parallel  under 
the  assumption  that  the  shape  parameter  is  independent  of  stress 
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The  purpose  of  the  test  discussed  in  [54]  was  to  determine  if  the  design 
requirement  of  a  20,000  hour  median  life  at  the  design  temperature  (180°C)  could 
be  met.  To  make  this  determination,  the  median  life  at  each  temperature  is  found  by 
noting  where  the  black  vertical  lines  in  Figures  16  and  17,  denoting  the  median 
failure  probability,  crosses  each  fit  line.  These  median  values  are  used  to  construct  a 
relationship  plot  (Figure  18)  showing  how  the  median  life  changes  with  increasing 
temperature.  Because  temperature  is  the  accelerating  variable,  the  Arrhenius 
reaction  rate  [45]  equation  was  used  to  model  the  life-temperature  relationship. 

The  Arrhenius  life-stress  model  is  given  by 

L(T)  —  C  expCAT-1)  (40) 

where  L(T )  represents  a  quantifiable  life  measure  such  as  characteristic  life  or 
median  life,  A  and  C  are  model  parameters  to  be  determined,  and  T  is  the  absolute 
temperature  in  degrees  Kelvin.  The  parameter  C  is  related  to  the  specimen 
geometry  and  testing  methodology,  while  A  is  a  function  of  the  activation  energy 
required  for  a  reaction  to  proceed  [54].  To  determine  if  the  Arrhenius  model 
accurately  describes  the  life-temperature  relationship  exhibited  by  the  data,  a 
process  similar  to  what  was  described  for  fitting  candidate  distributions  to  the  raw 
data  at  each  temperature  can  be  employed.  The  Arrhenius  model  is  first  "linearized” 
to  obtain 

In  [LCD]  =  ACT"1)  +  ln(C).  (41) 

If  the  Arrhenius  model  is  correct,  the  ln[L(T)]  values  plotted  against  the  reciprocal 
of  the  absolute  temperature  on  standard  plotting  paper  will  approximately  fall  on  a 
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line  with  slope  A  and  intercept  ln(C)  for  each  set  of  temperature  data.  Upon 
analyzing  the  data,  the  Arrhenius  parameter  values  for  the  plot  in  Figure  18  were 
estimated  to  be  A  —  6640.1  and  C  —  4.055  ■  10-3. 


Figure  18  -  Arrhenius  relationship  plot  of  the  example  data 

Alternatively,  the  data  may  be  presented  with  respect  to  temperature,  rather  that 
reciprocal  temperature  as  shown  in  Figure  19.  In  this  figure  the  median  times  to 
failure  (solid  line)  are  plotted  along  with  the  0.1  and  0.9  quantile  times  to  failure 
against  the  temperature  in  degrees  Celsius.  Extrapolating  the  solid  line  in  this  figure 
back  to  the  use  temperature  of  180°C  clearly  shows  that  the  current  median  life  of 
9000  hours  falls  well  short  of  the  20,000  hour  design  requirement. 

In  the  next  section  we  revisit  this  example,  utilizing  the  proposed  approach  in 
which  corrective  action  is  implemented  after  testing  is  performed  at  each 
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temperature  level.  The  impact  of  implementing  corrective  action  is  investigated 
assuming  the  times  to  failure  observed  at  220 °C  in  this  second  iteration  are 
identical  to  the  results  discussed  above.  A  benefit  of  the  proposed  approach  is  that 
the  number  of  required  samples  can  be  reduced  since  corrective  action  will  allow 
the  samples  to  be  reused  at  each  temperature  level. 


Figure  19  -  Adjusted  Arrhenius  relationship  plot  of  example  data  versus 

temperature 

Single  Failure  Mode  Case  -  With  Corrective  Action 

After  a  test,  each  failure  is  investigated  to  determine  its  root  cause  and  an 
appropriate  corrective  action  strategy  is  then  designed  to  eliminate  the  failure  in 
subsequent  tests.  However,  corrective  action  rarely  eliminates  a  failure  mode 
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completely.  It  is  common  in  the  reliability  growth  literature  [9]  to  characterize  the 
impact  of  a  corrective  action  by  modeling  the  reduction  in  a  failure  mode’s  hazard 
rate  by  a  proportional  amount  p  G  (0,1)  known  as  the  fix  effectiveness  factor  (FEF). 
The  FEF  varies  for  each  failure  mode,  but  [9,  p.  29]  reports  that  the  typical  range  of 
values  for  many  government  and  industry  systems  is  (0.55, 0.85)  with  a  typical 
value  of  0.70.  The  FEF  value  for  a  failure  mode  is  often  subjectively  estimated  after 
an  extensive  scoring  process  in  which  design  engineers  and  failure  analysis  experts 
consider  various  corrective  action  alternatives.  As  a  result  of  implementing 
corrective  action  the  system  reliability  improvement  can  be  expressed  as 


h'it) 

h(t) 


1  ~P 


(42) 


where  h'(t)  and  h(t)  denote  the  hazard  rate  functions  of  the  corrected  and 
uncorrected  systems,  respectively.  For  failure  modes  with  Weibull  distributed 
times  to  occurrence,  the  reduction  in  hazard  rate  resulting  from  corrective  action 
with  fix  effectiveness  factor  p  can  be  represented  as 

h\t )  fit?-1  ,a,P 


1  ~P  = 


h(t )  ( a 


—  ELI]  =(—Y 

')£  y  (a)^  J  Va'7 


(43) 


Substituting  the  Arrhenius  model  for  the  Weibull  scale  parameter  in  (47)  results  in 


,  P  P 

1  —  p  —  (CeA /CeA  )t  =  [exp(A  —  A')]t. 


(44) 


The  simplification  in  (44)  is  made  under  the  assumption  that  the  parameter  C  is 
unchanged  with  corrective  action.  Thus,  the  reduction  in  failure  rate  resulting  from 
a  corrective  action  with  fix  effectiveness  p  can  be  attributed  to  an  increase  in  the 
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activation  energy  required  to  initiate  the  failure  process,  where  the  increased 
activation  energy  is  expressed  as 

A'  —  A  —  (^)  ln[l  —  p\.  (45) 

However,  (45)  implies  that  the  value  of  A!  resulting  from  corrective  action  depends 
on  the  temperature  at  which  the  failure  was  discovered.  While  increased 
temperatures  allow  failures  to  be  observed  more  quickly,  the  strategy  to  remove  the 
failures  is  independent  of  temperature,  thus  there  exists  a  unique  value  of  A'  for 
each  corrective  action.  In  traditional  reliability  growth  models,  the  estimated  value 
of  p  is  based  on  an  assumed  mission  profile  of  the  fielded  system  [9,  p.  2]. 
Introducing  accelerated  stresses  to  a  reliability  growth  test  does  not  effect  this 
interpretation  of  p  and  the  improved  activation  energy  is  expressed  as 

A'  —  A  —  ln[l  -  p\.  (46) 

Table  13  displays  the  projected  reliability  measures  at  the  use  stress  of  180°C  that 
result  from  implementing  three  levels  of  corrective  action  after  the  first  test 
sequence. 


Table  13  -  Projected  reliability  measures  at  180°C  after  corrective  action 


FEF  (p) 

a ' 

t' 

Lmed 

A'  —  A 

0.85 

12855.95 

12094.14 

143.28 

0.70 

11453.35 

10774.65 

90.93 

0.55 

10704.93 

10070.59 

60.31 

Figure  20  displays  the  Arrhenius  plots  corresponding  to  the  projected  reliability 
measures  in  Table  13  along  with  the  original  plot  of  the  uncorrected  system.  Each 
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plot  in  the  figure  has  common  intercept  ln[C]  and  slope  increased  by  A'  —  A  over  the 
uncorrected  system.  The  vertical  line  in  the  figure  denotes  the  use  temperature 
180 °C.  The  ordinate  of  this  line  is  marked  to  indicate  the  20,000  hour  requirement. 


It  can  be  seen  in  Figure  20  that  none  of  the  projected  Arrhenius  plots  corresponding 
to  the  three  FEF  values  will  enable  the  system  to  meet  the  20,000  hour  requirement. 
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Figure  20  -  Updated  Arrhenius  plot  of  example  data  after  corrective  action 


This  same  conclusion  could  also  be  drawn  by  observing  the  projected  failure 
distributions  associated  with  each  of  the  three  levels  of  fix  effectiveness  p  as  shown 
in  Figure  21.  The  projected  distribution  function  WEIB(a' /?)  resulting  from 
corrective  action  is  expressed  as 
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P(X  <  x;  T)  =  1  —  exp 


C  exp 


In 


(1  ~P) 


1use 

p 


+  A\T 
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where  T  represents  any  temperature  of  interest.  Figure  21  presents  the  density  and 
reliability  functions  for  the  uncorrected  system  at  180°C  along  with  the  projected 
density  and  reliability  functions  after  corrective  action.  The  vertical  lines  spanning 
both  plots  in  the  figure  represent  the  median  values  for  each  distribution.  Again,  it 
is  clear  that  a  single  corrective  action  is  insufficient  to  achieve  the  requirement  as 
the  maximum  projected  median  life  is  just  over  12,000  hours. 


Figure  21  -  Joint  plot  of  density  and  reliability  functions  at  180°C  with 

corrective  action 
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Assuming  corrective  action  had  been  implemented  to  reduce  the  intensity  of  the 


failure  mode  discovered  during  the  220 °C  test,  the  improved  system  would  be  less 
sensitive  to  elevated  temperatures.  Quickly  discovering  additional  failures  in  the 
improved  system  would  therefore  require  the  application  of  a  more  severe 
environment  than  was  used  to  identify  the  failure  mode  initially.  For  systems  with  a 
single  failure  mode,  the  fix  effectiveness  of  a  second  corrective  action  may  be  far 
lower  than  that  of  the  initial  corrective  action.  Assuming  k  corrective  actions  are 
feasible  system  reliability  can  be  projected  recursively  and  the  Arrhenius  line  slope 
after  the  kth  corrective  action  becomes 


K 


A(fe)  =  ^  (in  (l  -p(/c))  k 

k= 1  ' 


+  A 


(48) 


and  the  projected  distribution  function  of  system  life  at  Tuse  is  expressed  as 

x 


/ 


P(X  <  x;  T)  =  1  —  exp 


V[ 


(  v( 

Tuse 

\  \  l 

C  exp 

(in 

J  +  dJr-i 

(49) 


) 


A  final  analysis  of  the  test  scenario  discussed  in  [54]  revealed  that  attaining  the 
20,000  hour  median  failure  time  at  the  use  stress  of  180°C  would  require  a  97% 
reduction  in  the  initial  hazard  rate.  A  reduction  of  this  magnitude  would  correspond 
to  implementing  three  corrective  actions  with  each  corrective  action  having  a  very 
high  fix  effectiveness  factor  of  0.85.  This  unlikely  event  confirms  the  additional 
commentary  provided  by  Nelson  [54,  p.  403]  where  further  testing  on  a  subsequent 
design  indicated  that  it  was  still  insufficient  to  attain  the  requirement  and  the 
project  was  abandoned.  Using  the  proposed  model  would  have  provided  an 
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indication  of  the  risk  in  attaining  the  required  reliability  while  utilizing  less  than  half 
of  the  specimens. 

Design  of  Accelerated  Reliability  Growth  Tests 

The  typical  objectives  in  designing  an  accelerated  reliability  growth  test  are  to 
minimize  the  total  test  time,  the  number  of  required  samples  and  the  prediction 
variance  of  the  estimate  of  system  life  in  the  usage  environment.  Because 
qualitative  testing  is  conducted  during  the  early  phases  of  product  development,  the 
number  of  prototypes  available  for  testing  is  already  severely  limited  due  to  the  cost 
of  production.  Thus,  identifying  a  best  qualitative  test  design  is  reduced  to  selecting 
the  stress  levels  at  which  the  accelerated  test  will  be  performed  and  how  the 
available  samples  will  be  allocated  to  those  stresses. 

To  model  lifetime  data  from  an  accelerated  test,  it  is  generally  assumed  that  the 
log  of  the  failure  times  observed  at  each  stress  level  follow  a  location-scale 
distribution 

P(Y  <  y)  =  <t>  (50) 

where  /r  and  o  are  the  true  but  unknown  location  and  scale  parameters  and  O  ( ■  ) 
denotes  the  standard  form  of  the  location-scale  distribution.  Further,  it  is  assumed 
that  the  value  of  the  scale  parameter  a  is  independent  of  stress  while  the  mean  log- 
failure  time  is  represented  as  a  linear  function  of  the  coded  stress  level  Xj  expressed 
as 

Ufa)  =  c0  +  c&j.  (51) 
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For  the  Arrhenius  relationship,  the  coded  stress  level  is  Xj  —  B /Tj  where  B  is 
Boltzmann’s  constant  and  7}  refers  to  the  temperature  in  an  absolute  scale  such  as 
Kelvin  or  Rankine.  The  variance  of  the  estimate  of  system  life  at  the  use-level  stress 
level  x0  as  result  of  observing  failures  at  stress  level  xb  i  —  1,2, ... ,  I  is  expressed  as 


Var[m(x0 )]  =  jl  +  (x0  —  x)2 


n 


m 


(52) 


where  m(  )  denotes  the  life  measure  of  interest,  such  as  the  mean  or  median,  and  n 
is  the  number  of  specimens  in  the  sample.  The  variance  in  (52)  is  minimized  by 
allocating  the  n  specimens  to  two  stress  levels,  a  high  stress  level  xH  and  a  low 
stress  level  xL  where  the  high  stress  level  is  defined  as  the  upper  limit  stress  level 
that  a  specimen  can  withstand  such  that  the  failure  modes  are  representative  of  the 
type  of  failures  expected  to  occur  in  the  usage  environment.  The  use  stress  level  x0 
and  high  stress  level  xH  establish  the  bounds  for  the  design  region  while  the  low 
stress  level  xL  is  chosen  to  minimize  the  variance  in  (52).  To  determine  the  optimal 
value  for  xL  in  a  general  design  region  Nelson  [54]  defines  the  extrapolation  factor  \ 


>  =  xH~xi 
XH  —  XL 


(53) 


At  the  use  stress  level  x0,  <f0  represents  the  ratio  of  the  allowable  design  range  xH  — 
x0  to  xH  —  x0  the  design  region  selected  for  the  accelerated  test.  It  can  be  shown 
that  the  variance  expression  in  (52)  can  be  simplified  as 

Ofo  -  vY 


Var[m(x0 )]  = 


1  + 


P(l-P) 


o 


n 


(54) 


when  two  test  stresses  are  chosen  that  lie  at  the  extremes  of  the  design  region,  that 
is  xH  and  xL,  where  p  denotes  the  proportion  of  specimens  allocated  to  the  low 
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stress  level  xL.  Further,  the  optimal  allocation  of  specimens  which  minimizes  the 
variance  in  (54)  for  a  general  design  region  is  denoted  by  p*  and  expressed  as 


V 


* 


2?o  -  1 ' 


(55) 


Substituting  p*  into  (54)  results  in  the  minimum  prediction  variance  of  the  estimate 


of  system  life  solely  as  a  function  of  the  extrapolation  factor 


Var*[m(x0 )]  =  1  +  4^0(^0  -  1)  I  —  I. 


In  the  extreme  case  of  xL  —  x0,  both  and  p*  equal  one  regardless  of  the  value  of  xH 
indicating  that  the  variance  is  minimized  when  all  n  tests  are  conducted  at  x0. 
However,  such  a  design  would  also  result  in  the  maximum  possible  test  duration.  As 
xL  increases,  the  optimal  allocation  of  samples  at  the  use  stress,  p*,  decreases  and 
reaches  a  minimum  value  of  0.5  as  <f0  ->  oo.  Clearly,  the  total  test  time  is  minimized 
by  raising  xL  and  xH  to  their  highest  possible  levels.  As  xH  and  x0  converge, 
however,  achieving  a  significant  reduction  in  test  time  requires  pushing  xH  closer  to 
the  destruct  limit.  Thus  care  must  be  taken  to  ensure  that  the  failures  observed  at 


the  high  stress  are  representative  of  the  type  of  failures  that  may  be  observed  in  the 
use  environment.  For  this  reason  the  two-stress  test  design  may  not  be  robust  to 
the  modeling  assumptions  and  the  addition  of  a  third  stress  level  may  be  beneficial. 
Figure  22  plots  the  log  variance  as  a  function  of  the  extrapolation  factor  for  the 
optimal  two-stress  test  design  along  with  the  two-  and  three-stress  level  designs 
with  an  equal  allocation  of  samples  at  each  stress  level. 
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When  the  allowable  design  region  xH  —  x0  is  large,  greater  freedom  exists  to  select 
xL  such  that  Var[m(x0 )]  can  be  minimized  while  simultaneously  reducing  the  total 
test  time.  Figure  23  shows  plots  of  the  total  test  time  against  the  xH/x0  ratio  for 
four  values  of  <f0.  For  each  plot,  x0  =  180°C,  as  was  the  case  in  the  example  test 
scenario  previously  described.  In  each  plot  the  low  stress  level  xL  increases  or 
decreases  commensurate  with  xH,  such  that  an  increase  in  xH  reduces  the  time  to 
failure  for  all  n  units  in  the  sample. 
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Figure  23  -  Comparison  of  the  reduction  in  total  test  time  for  the  two-stress 
equal  apportionment  design  the  two-stress  optimal  design  and  the  three- 
stress  equal  apportionment  design  for  various  levels  of  extrapolation  factor 


Identifying  a  best  design  for  a  given  test  scenario  with  inputs  x0  and  xH  could  be 
determined  through  multicriteria  optimization  where  the  objective  is  to 
simultaneously  minimize  Var[m(x0 )]  and  the  total  test  time  subject  to  the 
constraint  that  the  proportions  of  samples  tested  at  each  stress  level  sum  to  one. 
Figures  22  and  23,  indicate  that  the  three  stress  level  design  is  optimal  with  respect 
to  minimizing  the  total  test  time  but  sub-optimal  for  minimizing  the  variance 
compared  to  the  two-stress  level  designs  for  all  values  of  ^0.  Alternatively,  an 
optimal  design  minimizes  the  asymptotic  variance,  but  maximizes  the  total  time  on 
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test.  Therefore,  the  two-stress  equal  apportionment  design  is  recommended  as  a 
tradeoff  between  the  two-stress  optimal  design  and  the  three-stress  equal 
apportionment  design. 

Multiple  Failure  Modes  Case  -  With  Corrective  Action 

For  systems  containing  multiple  independent  failure  modes,  estimating  system 
reliability  from  accelerated  test  data  requires  sufficient  samples  to  observe  each 
failure  mode  at  multiple  stress  levels  [54].  Separate  models,  similar  to  those 
discussed  for  the  single  failure  mode  case,  may  then  be  developed  for  each  distinct 
failure  mode  and  the  exact  reliability  for  the  overall  system  can  then  be  determined 
as 

I 

R(t)  =  n^(0.  ;  =  1,2, (57) 

7  =  1 

In  qualitative  testing,  however,  it  is  rare  that  a  sufficient  number  of  prototypes  can 
be  made  available  to  account  for  all  of  the  failure  modes  that  have  yet  to  be 
discovered  in  early  system  testing.  For  tests  conducted  with  limited  sample  sizes, 
an  estimate  of  system  reliability  may  be  obtained  by  combining  the  observations  of 
multiple  failure  modes  to  form  a  single  distribution  plot  as  was  shown  above  for 
systems  with  only  one  failure  mode.  In  this  scenario  a  possibly  unknown  number  of 
independent  flaws,  denoted  by  J,  compete  to  be  the  cause  of  system  failure  (Figure 
24).  The  observed  lifetime  for  prototype  i  is  therefore  represented  as  the  minimum 
occurrence  time  among  the  J  modes  in  the  system. 
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Figure  24  -  Serial  arrangement  of  flaws  within  the  prototypes  subjected  to 
qualitative  accelerated  reliability  test.  The  arrangement  demonstrates  the 
competing  risk  assumption  in  the  model  where  the  time  to  failure  for 
prototype  i  is  the  minimum  activation  time  among  the  /  flaws. 

With  respect  to  flaw  j  E  {1, ... ,/}  the  probability  of  survival  up  to  time  t  is 
represented  by  the  random  variable  7}  with  sub-survivor  function  S(j,  t)  = 
Pr(Caitse  =  j,Time  >  t).  The  sub-survivor  function  is  not  a  proper  reliability 
function  in  that  S (J,  t)  =£  P(T  >  t\C  —  j ).  Instead  the  true  reliability  function 
specific  to  mode  j  is  expressed  as 


R(j,  0 


s(j,  t ) 

9; 


(58) 


where  qj  —  P(C  =  j )  =  F(j,  oo)  =  R(j,  0)  subject  to  the  constraint  2y=1  qj  —  1. 

The  form  of  R (J,  t)  in  (58)  may  be  unknown  prior  to  testing  and  is  often  assumed 
based  on  the  testing  scenario  or  prior  knowledge.  McLean  [80]  assumed 
exponentially  distributed  occurrence  times  for  each  failure  mode  on  the  basis  of 
mathematical  simplicity.  A  more  general  assumption,  suggested  for  the  proposed 
model,  is  R(j,  t)  ~  WEIB  ( ay,/?; ).  For  such  cases  where  J  >  2,  the  observed  failure 
time  Xt  —  min[7)]  follows  a  poly-Weibull  [85]  distribution  with  vector-valued 
parameters  a  and,  /?  and  density  function  expressed  as 
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(59) 


The  poly-Weibull  distribution  arises  naturally  in  scenarios  of  competing  risks  as  it 
describes  the  minimum  of  several  independent  random  variables  when  each  follows 
a  distinct  Weibull  law.  When  only  two  flaws  exist  (J  =  2),  the  distribution  is  known 
as  the  bi-Weibull,  and  J  —  3  model  is  naturally  called  the  tri-Weibull  distribution. 

An  advantage  of  the  poly-Weibull  distribution  lies  in  its  capacity  for  modeling 
not  only  increasing,  constant  and  decreasing  hazard  functions  but  also  non¬ 
monotone  hazard  functions  [95]  such  as  the  bathtub  curve  (Figure  25).  This 
property  is  important  as  non-monotone  hazard  functions  are  common  in  practice 
where  a  system  may  undergo  an  initial  "burn-in"  prior  to  periods  of  useful  life  and 
eventual  wearout. 


6=200,  50,  1.0,10 
6  =  100,250,  0  5,4.0 
6=  75,150,  2.5,5  0 
6  =  140,  40,10.0,1.2 
6=  36,160,  6.5,0  6 


o.oo- 


0  50  100  150  200  250  300 

t 

Figure  25  -  Example  poly-Weibull  hazard  functions  with  various  parameter 

values 
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If  the  shape  parameters  f31, ... ,  are  not  sufficiently  distinct  the  poly-Weibull 
may  over-fit  the  data  [27,  28],  For  the  extreme  case  where  /?1( ...  ,  /?y  =  /?,  the  poly- 

Weibull  converges  to  a  Weibull  distribution  with  density  /(t;  A,  /?)  =  A 
and  scale  parameter 


Our  interest  is  not  in  estimating  ccj,  /?;  for  all  j  G  {1, ... ,/},  as  the  data  produced  in 
qualitative  testing  is  insufficient  for  this.  Moreover,  the  poly-Weibull  CDF  cannot  be 
linearized,  as  was  done  for  the  Weibull  distribution  previously,  due  to  the 
summation  inside  the  logarithm  operator.  Thus,  the  probability  plotting  procedures 
discussed  above  for  the  single  failure  mode  case  cannot  be  employed  for  the  bi  and 
tri-Weibull  distributions.  Instead,  we  demonstrate  through  simulation  that  for 
scenarios  specific  to  qualitative  accelerated  life  testing  the  performance  of  Weibull 
distribution  is  sufficient  for  fitting  data  sets  in  which  the  observations  represent  the 
minimum  time  among  several  competing  failure  modes  each  having  Weibull 
distributed  occurrence  times. 

Monte  Carlo  Simulation  Procedure 

Prior  to  a  qualitative  reliability  test,  uncertainty  exists  regarding  the  number  of 
latent  failure  modes  embedded  within  each  prototype  and  the  expected  time  to 
occurrence  for  each  mode.  To  investigate  the  validity  of  the  Weibull  distribution  in 
the  multi-failure  mode  case,  a  Monte  Carlo  simulation  study  was  conducted.  The 
simulation  utilized  a  5  x  4  x  3  full  factorial  design  with  three  factors:  system 
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complexity,  sample  size,  and  initial  system  quality  as  depicted  in  Figure  26.  System 
complexity  is  represented  by  the  total  number  of  latent  failure  modes  competing  to 
be  the  cause  of  failure  in  each  identical  system  under  test.  For  this  study,  five  level 
levels  of  system  complexity  (J  —  5, 10, 20, 40, 50)  were  included.  Next,  four  levels  of 
sample  size  (n  =  10, 20, 40, 100)  were  utilized  to  illustrate  how  the  performance  of 
the  Weibull  distribution  compares  to  the  bi-  and  tri-Weibull  as  the  number  of  failure 
observations  increases.  Finally,  the  initial  system  quality  is  a  categorical  variable 
with  three  levels  (low,  middle,  high). 

We  suggest  that  products  with  high,  middle  and  low  levels  of  initial  quality  can 
have  an  identical  number  of  embedded  flaws  caused  by  poor  design  or 
manufacturing  defects.  However,  in  products  with  low  initial  quality  a  greater 
degree  of  variability  exists  in  the  severity  of  each  flaw  and  the  mean  time  until  a 
flaw  activates  may  be  far  lower  in  low  quality  products  than  for  a  similar  flaw  in 
products  with  higher  levels  of  quality. 
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Figure  26  -  Graphical  depiction  of  process  to  simulate  observed  minimum 
failure  times  for  systems  with  varied  numbers  of  competing  failure  modes  and 

levels  of  initial  design  quality 


To  simulate  the  time  to  occurrence  distribution  for  failure  mode  j  —  1, ...  ,J 
embedded  within  each  of  the  n  prototypes  in  the  sample  set,  a/  x  2  matrix  of  mode- 
specific  Weibull  parameter  values  was  created  for  each  sample.  The  Monte  Carlo 
approach  used  for  creating  these  failure  mode-specific  parameters  arises  from 
noting  that  for  random  variable  7}  ~  Weib(a.j,  fy'),  the  coefficient  of  variation  (Col/) 

7z 


CoV[T,]  = 


Var[Tj] 

v2 

r(i  +  2/ft)  1 

.  dtf . 

[r(i  + 1  /ft)2  J 

(61) 


is  solely  a  function  of  /?;  and  is  proportional  to  l//?;  as  shown  in  Figure  27.  Jin  et  al. 
[102]  notes  that  for  many  electronic  systems  the  CoV  tends  to  vary  between  0.05 
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and  5.50.  Thus,  an  approach  for  obtaining  the  mode-specific  Weibull  shape 
parameter  is  to  sample  a  CW[7)]  value  from  a  distribution,  substitute  this  value  into 


(61),  and  solve  for  jS ). 


0  5  10  15 

Pi 

Figure  27  -  Relationship  between  the  coefficient  of  variation  and  shape 
parameter  Bj  for  Weibull  distributed  random  variables 


The  support  region  for  this  sampling  distribution  should  mirror  the  range  defined 
by  Jin  et  al.  [102],  so  for  this  simulation  study  a  four-parameter  beta  distribution 
[103]  was  used.  IfX  ~  BETA(A,  6),  x  E  [0,1]  with  shape  parameters  A  and  6,  the 
four-parameter  beta  results  from  the  transformation  Y  —  X(d  —  c)  +  c  where  d  and 
c  are  the  upper  and  lower  limits  of  the  desired  support  for  the  transformed  variable 
Y.  Thus,  Y  ~  BETA{A,  9,  c,  d ),  y  E  [ c ,  d]  with  density  function  expressed  as 

1  (Y  —  c)A-1(d  —  T)0-1 


f(Y;A,  9,  c,  d)  = 


B(A,0)  (d-c) 


_ /'A  A+0 — 1 


(62) 
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The  combination  of  values  for  the  shape  parameters  A  and  0  can  then  be  used  to 
reflect  the  levels  of  design  quality  that  a  system  may  possess  at  the  start  of  a  test.  As 
illustrated  in  Figure  28,  the  combination  of  high  variability  Var  [Tj]  and  low  mean 
occurrence  times  E[Tj\  for  failure  modes  in  low  quality  systems,  imply  that  the 
likelihood  of  CW[7)]  values  at  the  upper  end  of  the  range  defined  by  [102]  is  greater 
than  in  systems  with  higher  quality  levels. 


0  1  2  3  4  5 

CoV 

Figure  28  -  Plots  illustrating  the  distribution  of  the  coefficient  of  variation 
values  for  multiple  levels  of  initial  design  quality 

In  this  simulation  study  the  sampling  distribution  for  CW[7)]  in  low  quality 
products  had  shape  parameters  (A,  0)  =  (1.25, 2.5)  to  reflect  an  increased  likelihood 
of  an  embedded  failure  mode  causing  infant  mortality.  To  model  products  with 
middle  and  high  levels  of  initial  quality,  0  =  5  and  10,  respectively,  were  used. 

These  parameter  values  result  in  a  sampling  distribution  for  CW[7)]  that  is  more 
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positively  skewed  than  that  of  low  quality  systems  indicating  a  reduced  likelihood  of 
infant  mortality  failures. 

The  Weibull  shape  parameter  obtained  from  the  sampling  procedure  outlined  in 
the  preceding  paragraphs  may  then  be  used  to  find  the  Weibull  scale  parameter  for 
failure  mode  j  by  substituting  /?;  into  the  expression 

OLj  =  E[Tj]  x  T(l  +  l/Pj)'1.  (63) 

As  previously  noted,  higher  coefficients  of  variation  imply  lower  mean 
occurrence  times,  therefore  the  sampling  distribution  for  E[Tj]  used  in  this  study  is 
conditioned  on  CW[7)].  For  this  simulation  E[Tj\  ~  BETA(y,  8, 3000, 0)  where  the 
support  region  (0,  3000)  represents  the  average  annual  usage  in  hours  of  an 


electronic  device.  The  corresponding  shape  parameters  y  and  8  are  determined  as 


(  CoVmax  -  CoV\T,]\  r  , 
y  =  1  +  rJT  ~  rJ  - 


CoVmax  —  CoVmin 


S  =  (i+C^M  e[1>4], 

V  CoVmax- CoVminJ 


Each  run  in  the  proposed  Monte  Carlo  simulation  can  be  parameterized  by  the 
vector  [ n,J ,  A,  6,  y,  5]  designating  the  sample  size,  system  complexity,  and  quality 
level  associated  with  a  given  test  scenario.  For  each  of  the  60  scenarios  the 
simulation  was  iterated  10,000  times.  During  each  iteration,  the  parameters 
cq  and  /?;,  j  =  1, ...  ,J  were  generated  and  the  J  Weibull  distributions  sampled  to 
provide  possible  failure  times.  From  these  J  values  the  minimum  is  chosen  as  the 
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observed  failure  time  for  prototype  i.  This  procedure  was  then  repeated  for  each  of 
the  remaining  prototypes  in  the  sample. 

Using  the  n  failure  observations  generated  in  each  sample,  the  relative 
performance  of  the  Weibull  distribution  was  compared  to  the  bi  and  tri-Weibull 
using  the  corrected  Akaike  Information  Criterion  ( AICc )  [104]  expressed  as 


2  k(k  +  1) 

AICc  =  2k  -  2(£)  +  — ^ - -  (65) 

n  —  k  —  1 

where  k  denotes  the  number  of  parameters  in  each  candidate  model  and  L  is  the 
value  of  the  log-likelihood  function  evaluated  at  the  maximum  likelihood  parameter 
estimates.  The  MLE’s  for  the  bi  and  tri-Weibull  distributions,  aMLE  and  pMLE,  were 
obtained  by  solving  the  system  of  2V  nonlinear  equations 

n 


dL 

dav 


v  =  l,...,V 


dL 

Wv 


V  =  1 . V 


(66) 


where  hv(t{)  =  (3vt ^  1ccJv  and  /i(u)  =  Y,v=i  lav^v)-  Solving  this  system  of 

equations  cannot  be  accomplished  analytically  and  numerical  root-finding 
techniques  can  be  tedious  as  finding  a  solution  can  be  sensitive  on  the  starting 
values  for  the  parameters  in  each  equation.  We  found  it  simpler  to  obtain  accurate 
parameter  estimates  by  maximizing  the  log-likelihood  function  directly  using  a 
quasi-Newtonian  method  [98]  and  utilizing  (66)  to  verify  the  solution.  For  each 
simulated  data  set,  the  AICc  of  one  of  the  three  candidate  distributions  has  the 
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lowest  value  implying  that  the  model  corresponding  to  AICcmin  has  the  highest 
relative  likelihood  of  minimizing  the  information  loss  expressed  by 

Kjriildata )  =  exp[0 .5(AICcmin  —  AlCc{j\.  (67) 

Appendix  A  displays  histograms  of  these  likelihood  values  for  the  Weibull,  bi- 
Weibull  and  tri-Weibull  distributions  for  each  of  the  sixty  testing  scenarios. 

Next,  Anderson-Darling  [105]  goodness  of  fit  tests  were  performed  to  determine 
the  likelihood  that  each  simulated  data  set  could  have  been  drawn  from  the 
WEIB(aMLE,  pMLE)  distribution.  The  unmodified  Anderson-Darling  test  is 
performed  using  the  test  statistic 

n 

A2  =  V - -  [ln(F(Yj)  +  ln(l  -  F(Yn+1_{)}  -  n.  (68) 

Z— <  n 

i 

where  F(  )  denotes  the  CDF  of  the  distribution  of  interest.  Stephens  [106]  found  A2 
to  be  one  of  the  best  empirical  distribution  function  based  test  statistics  for 
detecting  departures  from  several  continuous  distributions,  including  the  Weibull. 
Because  the  scale  and  shape  parameters  are  estimated  from  the  data  generated  in 
each  iteration  of  the  Monte  Carlo  simulation,  the  modified  Anderson-Darling 
statistic  [107]  A*  —  MA2  was  used.  Stephens  [107,  p.  146]  identified  the  functional 
form  of  M  for  the  cases  where  one  or  both  of  the  Weibull  parameters  are  estimated 
from  the  data.  Substituting  the  Weibull  distribution  function  and  M  —  1  +  yflrT1 
into  (72)  results  in  the  Anderson-Darling  test  statistic  used  for  this  study 
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For  each  iteration,  the  observed  significance  level  (p-value)  of  the  Anderson-Darling 
test  was  found  using  the  critical  value  tables  in  [107,  p.  146].  To  demonstrate  how 
the  Weibull  fit  the  simulated  data  generated  across  the  sixty  test  scenarios,  each 
iteration  of  the  model  was  considered  a  Bernoulli  trial  where  success  was  defined  as 
a  test  with  a  p-value  equal  to  or  greater  than  0.05.  Figure  29  displays  the  results  of 
the  analysis.  The  figure  shows  that  the  Anderson-Darling  test  fails  to  reject  the 
Weibull  for  all  regions  of  the  design  space  with  the  exception  of  the  low-quality, 
high-complexity  and  low-complexity,  high-quality  factor  level  combinations. 


10  Systems  - 20  Systems  - 40  Systems  - 100  Systems 


Figure  29  -  Results  of  the  Anderson-Darling  goodness  of  fit  tests  showing  the 
fraction  of  tests  with  significance  levels  greater  than  0.05  for  each 
combination  of  quality  level  sample  size  and  system  complexity 
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In  low-quality,  high-complexity  systems,  many  design  defects  or  other  potential 
causes  of  infant  mortality  failure  are  assumed  to  exist.  Thus,  the  simulated  failure 
observations  are  concentrated  in  a  time  interval  where  failure  occurs  almost 
instantaneously.  Fitting  data  sets  such  as  these  with  a  Weibull  distribution  implies  a 
shape  parameter  /?  «  1.  Observing  the  AICc  histograms  in  Appendix  A  for  systems 
with  these  factor  level  combinations  indicates  that  in  a  large  majority  of  the 
simulated  data  sets  the  Weibull  is  preferred  model.  Combining  the  results  displayed 
in  these  histograms  with  those  shown  in  Figure  29  suggests  that  the  low  quality 
level  factor  setting  may  be  too  low  to  represent  a  qualitative  accelerated  life  test. 

For  low-complexity,  high-quality  systems,  histograms  of  the  observations  indicate 
that  multimodal  density  functions  result  in  a  significant  proportion  of  the  generated 
data  sets.  It  is  in  this  design  region  that  non-monotone  hazard  rate  functions  such 
as  the  "bathtub”  curve  are  observed.  In  practice,  qualitative  accelerated  life  tests  are 
typically  performed  on  complex  electromechanical  devices  where  a  moderate  to 
high  level  of  initial  quality  is  expected. 

Observing  Figure  29  it  is  evident  that  for  systems  with  these  characteristics  the 
Weibull  distribution  is  adequate  to  model  an  overwhelming  majority  of  the 
simulated  data  sets.  This  result  implies  that  the  procedure  described  above  for 
projecting  system  reliability  after  corrective  action  in  the  single  failure  mode  case 
may  also  be  utilized  in  qualitative  testing  for  systems  with  multiple  mutually 
independent  failure  modes  where  each  follows  a  separate  Weibull  law. 
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Computing  the  system  FEF  when  multiple  competing  failure  modes  exist 
requires  taking  into  account  the  number  of  distinct  modes  exposed  during  the  test. 
The  Weibull  parameters  associated  with  each  failure  mode  are  unknown  as  is  the 
contribution  each  mode  makes  on  the  overall  failure  rate  of  the  system.  The  system 
failure  rate  function,  h(t)  =  /?t^  a~P ,  represents  the  sum  of  the  unknown  failure 
rates  of  each  of  the  embedded  modes  in  the  system  h{t )  =  YJj  hj(t ).  Assuming  each 
distinct  failure  mode  discovered  during  the  test  contributes  to  the  system  failure 
rate  equally,  the  system  FEF,  P,  can  be  expressed  as 


gm=l  Pm 
M 


(70) 


where  M  denotes  the  number  of  distinct  failure  modes  discovered  during  testing. 
For  the  example  scenario  used  to  describe  the  single  failure  mode  case,  temperature 
was  the  accelerating  variable  and  the  Arrhenius  model  was  used  to  describe 
relationship  between  median  life  and  temperature.  Expanding  this  functional 
relationship  for  the  multiple  failure  mode  case,  the  projected  system  reliability  after 
corrective  action,  with  respect  to  the  change  in  activation  energy  is  expressed  as 

i  Aisel 


A  =  In 


+  A. 


(71) 


And  the  projected  probability  density  function  of  system  life  at  Tuse  after  corrective 
action  is  expressed  as 


P{X<x-,T )  =  l-exp 


(  V 


X 

yC  exp 

Tuse 

+A) 

1  T_1 

|  1use 

(72) 
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Conclusions  and  Future  Work 


A  modified  accelerated  life  test  and  data  analysis  procedure  was  presented.  The 
proposed  approach  is  the  first  to  integrate  lifetime  acceleration  models  with 
reliability  growth  planning  to  estimate  reliability  improvement  after  implementing 
corrective  action.  The  procedure  was  assessed  by  first  considering  the  case  where 
the  systems  under  test  contained  a  single  dominant  failure  mode.  Utilizing  the 
results  from  a  published  accelerated  life  test  [18],  it  was  demonstrated  that  the 
likelihood  of  attaining  a  reliability  requirement  could  be  estimated  with  the 
proposed  procedure.  Further,  it  was  shown  that  by  using  the  proposed  procedure, 
an  indication  that  a  system  may  fail  to  reach  the  reliability  requirement  could  be 
obtained  by  testing  fifty  percent  less  samples  than  are  required  by  traditional 
quantitative  accelerated  life  test  methods.  This  result  suggests  that  a  significant 
reduction  in  the  time  and  cost  associated  with  system  testing  could  be  realized  as 
the  proposed  method  may  be  used  to  signal  a  need  to  reassess  a  product’s  design  or 
reallocate  testing  resources  to  avoid  unnecessary  maintenance  costs  or  an 
expensive  redesign  after  product  release. 

The  procedure  was  then  extended  to  the  case  where  each  system  under  test 
contained  multiple  independent  competing  failure  modes.  It  was  shown  that  if  the 
time  to  occurrence  distribution  for  each  respective  failure  mode  follows  a  distinct 
Weibull  law,  the  observed  system  failure  times  follow  a  poly-Weibull  distribution. 
However,  as  result  of  a  Monte  Carlo  simulation  it  was  shown  that  under  certain 
testing  scenarios,  specific  to  qualitative  accelerated  life  testing,  the  failure  data 
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generated  by  systems  with  multiple  independent  competing  failure  modes  can  be 
modeled  with  a  Weibull  distribution.  Thus,  the  proposed  model  utilizes  the  Weibull 
distribution  to  estimate  system  reliability  after  one  or  more  failure  modes  have 
been  discovered  and  removed. 
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VI.  Conclusions  and  Recommendations 


Dissertation  Summary 

Merging  reliability  growth  and  accelerated  life  testing  is  a  relatively  new  concept 
driven,  in  part,  by  the  complexity  and  widespread  use  of  modern  electronic  devices. 
In  commercial  environments,  where  several  product  types  may  be  in  development 
simultaneously,  little  time  may  exist  to  certify  that  the  reliability  of  a  particular 
product  meets  the  reliability  requirement.  Qualitative  accelerated  test  methods  are 
therefore  often  utilized  to  ensure  a  product’s  readiness  for  market  introduction  and 
to  avoid  costly  redesigns  caused  by  critical  failures  that  may  not  be  discovered  in 
traditional  reliability  testing.  While  published  case-studies  attest  to  the 
effectiveness  of  qualitative  life  testing  [15, 16, 14]  in  improving  product  reliabilty, 
the  capability  to  translate  the  limited  failure  data  into  a  meaningful  measure  of 
reliability  improvement  does  not  exist.  Thus  a  goal  of  this  dissertation  research  was 
to  bridge  the  gap  between  quantitative  and  qualitative  accelerated  reliability  test 
methods. 

In  Chapter  III,  several  commonly  used  accelerated  failure  data  analysis  methods 
were  highlighted  to  show  their  limitations  with  respect  to  estimating  system 
reliability  from  the  data  obtained  in  a  simple  qualitative  accelerated  stress  test.  It 
was  shown  that  these  methods  either  cannot  incorporate  design  changes  resulting 
from  corrective  action  or  over-simplify  the  reliability  growth  process  to  avoid 
estimating  the  parameters  of  the  life-stress  relationship.  It  was  demonstrated  that 
any  corrective  action  that  improves  a  product’s  reliability  changes,  by  definition,  the 
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way  the  product  responds  to  stress  and  therefore  also  alters  the  acceleration  factor. 
Further,  the  parameters  of  the  stress  life  relationship  must  be  estimated  for  each 
distinct  failure  mode  discovered  during  the  test.  Chapter  III  concluded  with  a 
discussion  on  several  potential  paths  for  future  research  with  respect  to  improving 
the  parameter  estimation  in  follow-on  testing. 

Because  accelerated  testing  models  correspond  to  individual  failure  modes,  the 
proper  approach  for  modeling  reliability  growth  of  a  complex  system  in  an 
accelerated  test  scenario  requires  an  accounting  of  the  failure  rate  changes  for  each 
mode  separately  and  integrating  these  changes  into  a  single  measure  of  system 
reliability.  For  each  failure  event,  only  one  risk  can  be  attributed  as  its  cause,  where 
a  risk  may  be  defined  as  either  a  single  failure  mode  or  the  combination  of  multiple 
modes  acting  simultaneously.  It  was  shown  in  Chapter  IV  that  if  it  can  be  assumed 
that  the  time  to  occurrence  distribution  for  each  competing  mode  follows  a  Weibull 
distribution  with  mode  specific  parameters  a;  and  /?;,  j  —  1, ... ,/,  the  system 
lifetime  can  be  modeled  with  the  poly-Weibull  distribution  with  vector  valued 
parameters  a,  /?.  The  poly-Weibull  was  then  presented  as  an  alternative  to  the 
Weibull  distribution  for  fitting  data  sets  where  the  hazard  function  is  monotone, 
non-monotone  or  even  bathtub  shaped.  A  numerical  and  analytical  procedure  for 
obtaining  the  poly-Weibull  maximum  likelihood  parameter  estimates  and 
asymptotic  standard  errors  was  presented  which  allows  reliability  engineers  to 
draw  conclusions  about  system  lifetime  with  a  stated  degree  of  confidence. 
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It  was  also  shown  that  the  poly-Weibull  distribution  is  capable  of  fitting 
reference  data  sets  known  to  have  bathtub-shaped  hazard  rate  functions  better  than 
the  best  modified  Weibull  models  in  the  literature.  The  goodness  of  fits  for  two 
forms  of  the  poly-Weibull  distribution,  the  bi-Weibull  and  tri-Weibull,  were 
assessed  using  the  Akaike  information  criterion  and  Kolmogorov-Smirnov  test 
statistic.  The  results  showed  that  the  bi-Weibull  and  the  tri-Weibull  outperform 
other  modified  Weibull  distributions  with  respect  to  their  fit  of  data  for  which  the 
hazard  rate  function  is  bathtub  shaped. 

Lastly,  a  model  was  presented  in  Chapter  V  to  characterize  a  product’s  likelihood 
of  attaining  a  reliability  requirement  after  implementing  corrective  action  to  remove 
one  or  more  embedded  failure  modes.  It  was  that  shown  under  the  assumptions  of 
qualitative  accelerated  life  testing  the  Weibull  distribution  may  be  used  to  model 
accelerated  life  test  data  when  multiple  independent  failure  modes  are  observed. 

An  example  test  scenario  [18]  was  discussed  in  which  a  reliability  enhancement 
program  was  abandoned  after  multiple  tests  showed  that  the  new  design  would  not 
meet  the  required  median  life.  Our  results  indicated  that  had  our  proposed  model 
been  used  to  estimate  the  reliability  risk,  a  significant  cost  and  time  savings  could 
have  been  realized  by  using  less  than  half  of  the  samples  required  for  the  original 
test.  Further,  in  Chapter  V  it  was  shown  that  the  projection  model  can  be  used  to 
conduct  trade  off  analyses  as  a  basis  for  reviewing  the  reliability  requirements. 
These  are  important  contributions  as  the  model  allows  reliability  engineers  to 
design  better  systems  at  less  cost  and  in  less  time  than  would  be  possible  otherwise. 
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Conclusions  and  Future  Work 


Qualitative  accelerated  life  test  methods  were  initially  developed  with  the  intent 
of  discovering  and  removing  previously  unknown  failure  modes  from  early  system 
prototypes.  As  opposed  to  quantitative  accelerated  life  testing,  these  methods  are 
not  designed  with  life  estimation  or  prediction  in  mind.  In  a  few  of  these  tests,  such 
as  highly  accelerated  stress  testing,  extensive  testing  and  data  analysis  has  been 
presented  in  the  literature  to  derive  a  life-stress  relationship  for  specific 
combinations  of  stressors.  In  general,  however,  the  ability  to  devise  ever-harsher 
test  environments  to  discover  more  failure  modes  in  less  time  will  always  outpace 
the  rate  at  which  mathematically  sound  tools  can  be  generated  to  estimate  or 
predict  product  reliability. 

The  goal  of  this  dissertation  was  to  provide  first  step  toward  bridging  the  gap 
separating  quantitative  and  qualitative  accelerated  life  testing.  The  most  obvious 
next  steps  to  be  taken  in  further  research  to  improve  the  reliability  growth 
projection  model  involve  extending  the  test  scenario  to  multiple  stressors  and  to 
include  multiple  corrective  action  periods  to  remove  failure  modes  discovered  in 
subsequent  tests.  Both  of  these  paths  are  discussed  in  the  sections  below. 

Extending  the  Model  to  Multiple  Corrective  Action  Periods 

For  systems  with  multiple  failure  modes,  it  is  unlikely  that  a  single  test  with  a 
single  corrective  action  period  can  uncover  all  of  the  failure  modes  that  could  occur 
in  the  use  environment.  In  reality,  many  corrective  action  periods  may  be  required. 
Extending  the  projection  model  to  testing  scenarios  with  multiple  corrective  action 
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periods  requires  that  the  corrective  actions  to  remove  the  failure  modes  discovered 
during  the  preceding  tests  are  not  cost  prohibitive.  The  reliability  growth  attained 
through  corrective  action  in  an  accelerated  test  depends  on  the  likelihood  of  the 
failure  actually  occurring  in  operation  [54,  79, 17].  For  the  single  failure  mode  case, 
it  was  shown  in  Chapter  V  that  each  corrective  action  increases  the  activation 
energy  by  A'  —  A.  With  each  successive  improvement,  the  energy  required  to 
activate  a  failure  mode  approaches  the  energy  required  for  a  non-representative 
failure  to  occur.  Nelson  [54]  describes  an  example  in  which  a  costly  program  was 
initiated  to  remove  a  failure  mode  that,  it  was  later  determined,  would  have  never 
occurred  in  the  use  environment.  Thus,  the  end  result  of  the  program  was  a 
substantial  amount  of  wasted  resources  and  no  actual  reliability  improvement. 

One  possible  way  to  address  this  is  by  introducing  an  adjusted  FEF  of  the  form 

p*=p-d(S )  (73) 

where  p  is  the  estimated  FEF  value  for  a  given  failure  mode  and  d(S)  G  (0,1)  is  a 
multiplicative  factor  that  scales  p  based  on  the  level  of  stress  S  at  which  the  test  is 
performed.  When  the  stress  applied  to  the  system  is  at  or  below  the  field-use  stress, 
S  <  Suse,  d(S)  is  defined  to  be  equal  to  unity.  Conversely,  when  the  applied  stress  is 
above  some  limit  stress  S  >  Siimit, where  failure  occurs  almost  immediately  due  to 
an  unrecoverable  failure,  d(S)  =  0.  Examples  of  unrecoverable  failures  include  the 
melting  or  phase  change  of  a  material  as  result  of  elevated  temperatures, 
catastrophic  mechanical  failure  due  to  elevated  vibrational  stresses  and  dielectric 
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breakdown  due  to  voltage  overstress.  For  Suse  <  S  <  Siimit,  a  possible  form  of  d(S) 
is  a  monotone  decreasing  sigmoid  function  such  as 


d(S)  =  1 


1  +  exp  [-  L - w_  )  (5  -  (5ttsc 

L  ^ limit  ^ use  '  \  ^  ) 


(74) 


where  the  functional  form  of  (74)  is  derived  by  modifying  the  logistic  function.  The 
parameter  w  may  be  used  to  represent  how  quickly  d(S)  approaches  zero  for  a 
given  system  and  testing  scenario.  Figure  30  illustrates  the  shape  of  d(S)  for 
various  values  of  Siimit  and  for  w  —  10.  The  true  form  of  d(S)  is  unknown,  however 
the  figure  illustrates  that  accelerated  reliability  growth  test  implies  a  trade-off 
between  the  rate  at  which  failures  occur  and  the  level  of  improvement  made.  By 
introducing  an  adjusted  FEF  value  to  the  model  developed  in  Chapter  V,  test 
planners  can  balance  the  level  of  stress  applied  with  the  expected  reliability 
improvement. 


Figure  30  -  Sigmoid  shape  of  d(S) 
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Extending  the  Model  to  Multiple  Stressor  Tests 

Extending  the  reliability  growth  projection  model  developed  in  Chapter  V  to 
address  stress  states  in  which  more  than  one  stressor  is  used  would  require 
estimating  the  parameters  of  the  life-stress  relationship  for  each  stress  separately 
along  with  any  interaction  terms  that  may  be  significant.  Even,  if  it  can  be  assumed 
that  no  interaction  exists  between  the  stressors  (a  very  restrictive  assumption)  each 
stressor  added  to  a  design  would  bring  two  additional  parameters  that  must  be 
estimated.  Estimating  these  parameters  would  require  multiple  observations  of 
each  failure  mode  from  tests  performed  at  distinct  design  points.  By  dividing  an 
already  small  test  sample  across  a  greater  number  of  design  points  would  likely 
result  in  parameter  estimates  of  the  life-stress  relationship  that  have  large 
variances.  This  would  lead  to  confidence  intervals  that  would  be  far  too  wide  to  be 
meaningful.  Unfortunately,  financial  realities  limit  the  number  of  units  that  can  be 
made  available  for  early  system  testing  and  therefore  constrain  the  decision  space 
as  to  feasible  testing  options.  It  is  our  hope  that  in  spite  of  these  constraints 
progress  can  be  made  to  further  bridge  the  gap  between  qualitative  and  quantitative 
accelerated  life  testing. 
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Appendix  A:  Poly-Weibull  Observed  Fisher  Information  Equations 

The  observed  Fisher  information  for  the  bi-Weibull  and  tri-Weibull  distributions 
are  expressed  respectively  as 


and 


0Bw  (a,/?)  =  — VVT£(a,  /?)  = 


0TW(a,p)  =  — VVT£(a,  /?)  = 


d2L 

m 

d2L 

9ft9ft 

d2L 

9ft9ft 

d2L 


da1dpi 

d2L 

da2dp1 

d2L 


L  da3dpi 


■  a2x 

d2L 

d2L 

d2L  - 

aft2 

d2L 

3ft3ft 

d2L 

dpida1 

d2L 

dP\da2 

d2L 

dp2dpl 

d2L 

dPl 

d2L 

dp2da1 

d2L 

dp2da2 

d2L 

da-^dP-L 

d2L 

da1dp2 

d2L 

da2 

d2L 

da1da2 

d2L 

_dcc2dp1 

da2dp2 

da2da1 

da2  \ 

d2L 

d2L 

d2L 

d2L 

d2L 

dpidp2 

d2L 

3ft  3ft 
d2L 

dpida1 

d2L 

9ft3a2 

d2L 

dPida3 

d2L 

3ft2 

d2L 

3ft3ft 

d2L 

dp2doc1 

d2L 

dp2da2 

d2L 

dp2da3 

d2L 

3ft  3ft 
d2L 

9ft2 

d2L 

9ft  9^ 
d2L 

dp3da2 

d2L 

dp3da3 

d2L 

3«i  3ft 
d2L 

3ai3ft 

d2L 

da2 

d2L 

da1da2 

d2L 

da1da3 

d2L 

da2dp2 

d2L 

da2dp3 

d2L 

da2da1 

d2L 

dal 

d2L 

da2da3 

d2L 

da3dp2 

da3dp3 

da3da1 

da3da2 

da2 

The  elements  of  these  matrices  are  listed  below,  where  j,  k  —  1,2,  3, ... 


dL 

Wj 


N 

^  hj(ti) 


(log [ti/aj)  +  Pj  ft 
hit,) 


d2L  ^  hj(ti)  ■  log (aj/ti)  (log(a//ti)  -  2ft  ft  ■  ft(tj)  -  (h;(ti)(log(ti/ay)  +  Pj  ft) 

Wd2 


d2£  _  ^  fty(ti)(log(ti/a;)  +  ft  1)  ■  ftk(ti)Oog(ti/ak)  +  ft 
dpjdpk  ~  4 

J  l 


a2x 

dPjdctj 


(-/i(ti)ft;(ti)(log(ti/a/  )  +  2ft  ft  +  hj(tt)2(log(tt/aj)  +  ft  ft) 


127 


dPjdak 


\  ’  /?.- 
/  hj(t{)hk{t{)  — 
L-i  1  CCj 


Pj  [Qog (ti/aj)  +  /?;  1)' 


JV 

3£  ^  Pj  tj  1 

dctj  l— i  cij  Pj  h(ti) 


d2L  ST' Pj  hj(ti)h(ti)(Pj  +  l)  ( Pjhj 


08  +  D 


dLh.r^A.Pk 


d*L  $  ■%**<.«) 


dajdak 


=  1- 


128 


Appendix  B:  Histograms  of  Corrected  Akaike  Information  Criterion  Values 
(AICc)  for  the  Weibull,  Bi-Weibull  and  Tri-Weibull  Distributions 

Using  the  n  observations  in  each  of  the  simulated  samples  described  in 
Chapter  V,  the  relative  performance  of  the  Weibull  distribution  was  compared  to  the 
bi-  and  tri-Weibull  using  the  Akaike  Information  Criterion  corrected  for  sample  size 
(AICc)  [104].  The  corrected  Akaike  Information  Criterion  is 

2  k(k  +  1) 

AICc  =  2k  —  2(L)  +  — ^ - -  (75) 

n-k-1 

where  k  denotes  the  number  of  parameters  in  each  candidate  model  and  L  is  the 
value  of  the  log-likelihood  function  evaluated  at  the  maximum  likelihood  parameter 
estimates.  For  each  simulated  data  set,  the  AICc  of  one  of  the  three  candidate 
distributions  will  have  the  lowest  value.  The  model  corresponding  to  AICcmin  has 
the  highest  relative  likelihood  of  minimizing  the  information  loss.  This  likelihood  is 
Kjriildata)  =  exp[0 .S(AICcmin  —  AICc{)]  G  [0,1].  (76) 

This  appendix  displays  histograms  of  the  likelihood  values  obtained  for  the 
Weibull,  bi-Weibull  and  tri-Weibull  distributions  for  each  of  the  sixty  testing 
scenarios.  In  each  figure  the  histogram  for  the  Weibull  distribution  is  presented  in 
the  top  plot  while  the  bi-Weibull  and  tri-Weibull  are  shown  in  the  middle  and 
bottom  plots,  respectively.  The  figures  are  organized,  first  according  to  the  number 
of  embedded  failure  modes,  then  by  the  number  of  systems  in  each  sample  and 
finally  by  the  quality  level  of  each  system. 
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Figure  31  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  10  systems,  quality  =  low) 
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Figure  32  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  20  systems,  quality  =  low) 
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Figure  33  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  40  systems,  quality  =  low) 
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Figure  34  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  100  systems,  quality  =  low) 
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Figure  35  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  10  systems,  quality  =mid) 
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Figure  36  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  20  systems,  quality  =  mid  ) 
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Figure  37  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  40  systems,  quality  =  mid  ) 
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Figure  38  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  100  systems,  quality  =mid) 
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Figure  39  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  10  systems,  quality  =high) 
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Figure  40  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  20  systems,  quality  =  high  ) 
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Figure  41  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  40  systems,  quality  =high) 
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Figure  42  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  5  modes,  n  =  100  systems,  quality  =high) 
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Figure  43  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  10  systems,  quality  =low) 
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Figure  44  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  20  systems,  quality  =low) 
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Figure  45  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  40  systems,  quality  =low) 
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Figure  46  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  100  systems,  quality  =low) 
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Figure  47  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  10  systems,  quality  =mid) 
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Figure  48  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  20  systems,  quality  =  mid  ) 
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Figure  49  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  40  systems,  quality  =mid) 
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Figure  50  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  100  systems,  quality  =mid) 
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Figure  51  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  10  systems,  quality  =high) 
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Figure  52  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  20  systems,  quality  =high) 
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Figure  53  -  Histogram  of  AlCc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  40  systems,  quality  =high) 
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Figure  54  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  10  modes,  n  =  100  systems,  quality  =  high  ) 
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Figure  55  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  10  systems,  quality  =low) 
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Figure  56  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  20  systems,  quality  =low) 
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Figure  57  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  40  systems,  quality  =low) 
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Figure  58  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  100  systems,  quality  =low) 
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Figure  59  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  10  systems,  quality  =mid) 
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Figure  60  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  20  systems,  quality  =  mid  ) 
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Figure  61  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  40  systems,  quality  =mid) 
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Figure  62  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  100  systems,  quality  =mid) 
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Figure  63  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  10  systems,  quality  =high) 
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Figure  64  -  Histogram  of  AlCc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  20  systems,  quality  =high) 
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Figure  65  -  Histogram  of  AlCc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  40  systems,  quality  =high) 
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Figure  66  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  20  modes,  n  =  100  systems,  quality  =high) 
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Figure  67  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  10  systems,  quality  =low) 


o 

c 

0 

=3 

cr 

0 


166 


Tri-Weibull  Bi-Weibull  Weibull 


frequency  frequency  frequency 


8000-1 


6000- 


2000- 


0j  - = — -  1  — - - — - - 

I - 1 - 1 - 1 - 1 - 1 

0.0  0.2  0.4  0.6  0.8  1.0 


Figure  68  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  20  systems,  quality  =low) 
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Figure  69  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  40  systems,  quality  =low) 
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Figure  70  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  100  systems,  quality  =low) 
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Figure  71  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  10  systems,  quality  =mid) 
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Figure  72  -  Histogram  of  AlCc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  20  systems,  quality  =  mid  ) 
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Figure  73  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  40  systems,  quality  =  mid  ) 
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Figure  74  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  100  systems,  quality  =mid) 


173 


Tri-Weibull  Bi-Weibull  Weibull 


8000-1 


o  6000- 

c 

0 

D 

4000- 


2000- 


0 

0.0  0.2  0.4  0.6  0.8  1.0 


10000-1 


8000- 

>. 

^  6000- 

0 
D 
CT 

0  4000- 

2000- 

0J - 

1 - 1 - 1 - 1 - 1 - 1 

0.0  0.2  0.4  0.6  0.8  1.0 

Figure  75  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  10  systems,  quality  =high) 
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Figure  76  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  20  systems,  quality  =high) 
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Figure  77  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  40  systems,  quality  =high) 
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Figure  78  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  40  modes,  n  =  100  systems,  quality  =high) 
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Figure  79  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  10  systems,  quality  =low) 
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Figure  80  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  20  systems,  quality  =low) 
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Figure  81  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  40  systems,  quality  =low) 
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Figure  82  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  100  systems,  quality  =low) 
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Figure  83  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  10  systems,  quality  =mid) 
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Figure  84  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  20  systems,  quality  =mid) 
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Figure  85  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  40  systems,  quality  =  mid  ) 
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Figure  86  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  100  systems,  quality  =mid) 
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Figure  87  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  10  systems,  quality  =high) 
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Figure  88  -  Histogram  of  AlCc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  20  systems,  quality  =high) 
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Figure  89  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  40  systems,  quality  =high) 
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Figure  90  -  Histogram  of  AICc  values  for  the  Weibull,  bi-Weibull  and  tri- 
Weibull  distributions  (/  =  50  modes,  n  =  100  systems,  quality  =  high  ) 
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Appendix  C:  Computer  Code  Used  to  Compare  the  Fit  of  Various  Distributions 
to  the  Aarset  and  Meeker  Datasets  as  Discussed  in  Chapter  IV 


#####  Meeker  Data  set 
Meekerc- 

sort(c(2  75, 13, 147, 22. 5, 23. 5, 181, 30, 65, 10, 173, 106, 2 12, 2, 261, 293, 88, 247, 28, 143, 80, 
245,266,300,300,300,300,300,300,300,300)) 

Meekerlc- 

sort(c(275, 13, 147, 22.5, 23.5, 181, 30, 65, 10, 173, 106, 212, 2, 261, 293, 88, 247, 28, 143, 80, 

245,266))#300,300,300,300,300,300,300,300)) 

Meeker2<-c(300,300,300,300,300,300,300,300) 


Censor<-c(l,  1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,1, 1,0, 0,0, 0,0, 0,0,0) 


eCDF<- 


sort(c(275, 13, 147, 22. 5, 23.5, 181, 30, 65, 10, 173, 106, 212, 2, 261, 293, 88, 247, 28, 143, 80, 
245,266,300)) 

RelFreqc- 

c(0.03333, .06667, 0.1, .13333, .16667, .2, .2333, .26667, .3, .3333, .36667, .4, .43333, .4666 
7, .5, .53333, .56667, .6, .63333, .66667, .7, .73333,1) 


#####  Aarset  Data  set 
Aarsetc- 

c(0.1, 0.2, 1,1, 1,1, 1,2, 3, 6, 7, 11, 12, 18, 18, 18, 18, 18, 21, 32, 36, 40, 45, 46, 47, 50, 55, 60, 63, 63, 
67,67,67,67,72,75,79,82,82,83,84,84,84,85,85,85,85,85,86,86) 
Aarset2<-c(Aarset+rnorm(50,  0,  0.001)) 

RelFreq<-c(rep(l/length(Aarset2),length(Aarset2)))*c(l:length(Aarset2)) 


Red<-c(0.1,  0.2, 

1,2,3,6,7,11,12,18,21,32,36,40,45,46,47,50,55,60,63,67,72,75,79,82,83,84,85,86) 

RedAA<- 

c(.02, .04, .14, .16, .18, .2, .22, .24, .26, .36, .38, .4, .42, .44, .46, .48, .5, .52, .54, .56, .6, .68, .7, .72, .74, 
.78, .8, .86, .96,1) 


#########  Compare  Fit  of  Poly-Weibull/New  Modified  Weibull  Reliability  Against 
Meeker  Data  ########### 


########  New  Modified  Weibull  Model  ######################### 
B0<-matrix(NA,  nrow=100,  ncol=5,  byrow=FALSE,);  colnames(B0)<- 
c("a","b","A","g","l") 

B0[,l]<-c(runif(100,  0.01,1))  #a 
B0[,2]<-c(runif(100,  0.01,1))  #b 
B0[,3]<-c(runif(100,  0.01,1))  #t 
B0[,4]<-c(runif(100,  0.01,1))  #g 
B0[,5]<-c(runif(100,  0.01,1))  #1 
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NMW<-function(x,y,d){ 

a<-x[l] 

b<-x[2] 

t<-x[3] 

g<-x[4] 

l<-x[5] 

n<-length(Meeker) 

F<-c(rep(NA,l)) 

F<-(sum(log(b*(g+l*y)*(yA(g-l))*exp(l*y)+a*t*yA(t-l))*d)-a*sum(yAt)- 

b*sum((yAg)*exp(l*y))) 

return(-F) 

} 

N<-matrix(NA,  nrow=100,  ncol=6,  byrow=FALSE,);  colnames(N)<- 
c("a","b";,t";,g","l", "Value") 

for  (i  in  1:100){ 

N[i,]=0 

N[i,l:5]<-optim(par=B0[i,],  fn=NMW,  y=c(Meeker),d=Censor,  method="L-BFGS-B", 
lower=c(0.0001, 10A-9,  0.0001,  0.0001,  0.0001),  upper=c(l,  1, 1, 1, 1))  $par 
N[i,6]<-NMW(N[i,l:5],c(Meeker),d=Censor) 

} 

NM<-which.min(N[,6]) 

####New  Modified  Weibull  PDF 
NMWPdf<-function(t){(N[NM,l]*N[NM,3]*tA(N[NM,3]- 
l)+(N[NM,2]*(N[NM,4]+N[NM,5]*t)*tA(N[NM,4]-l)*exp(N[NM,5]*t)))*(exp(- 
l*(N[NM,l]*tAN[NM,3]+N[NM,2]*tAN[NM,4]*exp(N[NM,5]*t))))} 

####New  Modified  Weibull  Reliability  function 
NMWeib<-function(t){exp(- 

l*(N[NM,l]*tAN[NM,3]+N[NM,2]*tAN[NM,4]*exp(N[NM,5]*t)))} 

####New  Modified  Weibull  Hazard  Function 

NMWHaz<-function(t){(N[NM,l]*N[NM,3]*tA(N[NM,3]- 

l)+(N[NM,2]*(N[NM,4]+N[NM,5]*t)*tA(N[NM,4]-l)*exp(N[NM,5]*t)))} 

####  New  Modified  Weibull  CDF 
NMcdf<-function(x,y){ 
a<-x[l] 
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b<-x[2] 

t<-x[3] 

g<-x[4] 

l<-x[5] 

F<-c(rep(NA,l)) 

F<-(l-exp((-a*yAt)-(b*yAg)*exp(l*y))) 

return(F) 


} 

NMtest<-NMcdf(x=c(N[NM,l],N[NM,2],N[NM,3],  N[NM,4],N[NM,5]),y=Meeker) 


###  Bi-Weibull  Model  ######################### 

pO<-matrix(NA,  nrow=100,  ncol=4,  byrow=FALSE,);  colnames(pO)<- 

c("betal","beta2","thetal","theta2") 

p0[,l]<-c(runif(100,  0.25,5 ))  #betal 

p0[,2]<-c(runif(100,  0.25,5 ))  #beta2 

p0[,3]<-c(runif(100,  min(Meeker),max(Meeker)))  #thetal 

p0[,4]<-c(runif(100,  min(Meeker),max(Meeker)))  #theta2 

GG<-function(x,y,d){ 

bl<-x[l] 

b2<-x[2] 

tl<-x[3] 

t2<-x[4] 

F<-c(rep(NA,l)) 

F<-sum(log(((bl*(yA(bl-l))*(tl)A(-bl))+(b2*yA(b2-l))*((t2)A(-b2))))*d- 

CC(y)/tl)A(bl)+(y/t2)A(b2))) 

return(-F) 

} 

P<-matrix(NA,  nrow=100,  ncol=5,  byrow=FALSE,);  colnames(P)<- 
c("betal","beta2","thetal","theta2",  "Value") 

for  (i  in  1:100){ 

P[iJ=0 
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P[i,l:4]<-optim(par=pO[i,],  fn=GG,  y=c(Meeker),d=c(Censor),  method="L-BFGS-B", 
lower=c(0.25,  0.25, 10, 10),  upper=c(90,  90,  400,  400))  $par 
P[i,5]<-GG(P[i,l:4],y=c(Meeker),d=c(Censor)) 

} . 

PM<-which.min(P[,5]) 

P[PMJ 

#SEB<-optim(par=P[PM,l:4],  fn=GG,  y=c(Meeker),d=l,  method="L-BFGS-B", 
lower=c(0.25,  0.25, 10, 10),  upper=c(90,  90, 100, 100),hessian=TRUE)  $hessian 
#seb<-sqrt(diag(solve(SEB))) 

####Bi-Weibull  Reliability  function 

BiWeib<-function(t){exp(-l*(((t/(P[PM,3]))AP[PM,l])+((t/(P[PM,4]))AP[PM,2])))} 

####Bi-Weibull  Hazard  function 

BiHaz<-function(t){(P[PM,l]*tA(P[PM,l]- 

1))/(P[PM,3]  AP[PM,l])+(P[PM,2]*tA(P[PM,2]-l))/(P[PM,4]  AP[PM,2])} 
####Bi-Weibull  pdf  function 
Bipdf<-function(t){((P[PM,l]*tA(P[PM,l]- 

1))/(P[PM,3]  AP[PM,l])+(P[PM,2]*tA(P[PM,2]-l))/(P[PM,4]  AP[PM,2]))*(exp(- 
l*(C(t/(P[PM,3]))AP[PM,l])+C(t/(P[PM,4]))AP[PM,2]))))} 

####Bi-Weibull  CDF 
BWcdf<-function(x,y){ 
bl<-x[l] 
b2<-x[2] 
tl<-x[3] 
t2<-x[4] 

F<-c(rep(NA,l)) 

F<-l-(exp(-l*(((y/tl)Abl)+((y/(t2))Ab2)))) 

return(F) 


} 

BWtest<-BWcdf(x=c(P[PM,l],P[PM,2],P[PM,3],P[PM,4]),y=Meeker) 


###########################  Tri-Weibull  Model 
######################### 

tO<-matrix(NA,  nrow=100,  ncol=6,  byrow=FALSE,);  colnames(tO)<- 
c("betal","beta2","beta3","thetal","theta2","theta3") 
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t0[,l]<-c(runif(100,  0.25,5 ))  #betal 
t0[,2]<-c(runif(100,  0.25,5 ))  #beta2 
t0[,3]<-c(runif(100,  0.25,5 ))  #beta3 
t0[,4]<-c(runif(100,  min(Meeker),max(Meeker)))  #thetal 
t0[,5]<-c(runif(100,  min(Meeker),max(Meeker)))  #theta2 
t0[,6]<-c(runif(100,  min(Meeker),max(Meeker)))  #theta3 


TT  <-function(x,y,d){ 
bl<-x[l] 
b2<-x[2] 
b3<-x[3] 
tl<-x[4] 
t2<-x[5] 
t3<-x[6] 

F<-c(rep(NA,l)) 

F<-sumnogffbl*(yAfbl-l))*ftl)A(-bl))+fb2*yAfb2-l))*(ft2)A(-b2))+fb3*yAfb3- 

l))*((t3)A(-b3)))*d-(((y)/tl)A(bl)+(y/t2)A(b2)+(y/t3)A(b3))) 

return(-F) 

} 

Tr<-matrix(NA,  nrow=100,  ncol=7,  byrow=FALSE,);  colnames(Tr)<- 
c("betal","beta2","beta3","thetal","theta2","theta3", "Value") 

for  (i  in  1:100){ 

Tr[i,]=0 

Tr[i,l:6]<-optim(par=t0[i,],  fn=TT,  y=c(Meeker),d=c(Censor),  method="L-BFGS-B", 
lower=c(25,  0.25,  0.25, 10, 10, 10),  upper=c(124,  90,  90, 1200,  400,  400))  $par 
Tr[i,7]<-TT(Tr[i,l:6],y=c(Meeker),d=c(Censor)) 

} 

TM<-which.min(Tr[,7]) 

Tr[TM,] 

#SET<-optim(par=Tr[TM,l:6],  fn=TT,  y=c(Meeker),d=l,  method="L-BFGS-B", 
lower=c(0.25,  0.25,  0.25, 10, 10, 10),  upper=c(110,  90,  90, 100, 125, 
100),hessian=TRUE)  $hessian 
#set<-sqrt(diag(solve(SET))) 

####Tri-Weibull  Reliability  function 
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TriWeib<-function(t){exp(- 

l*(((V(Tr[TM,4]))ATr[TM,l])+((t/(Tr[TM,5]))ATr[TM,2])+((t/(Tr[TM,6]))ATr[TM,3 

])))} 

####Tri-Weibull  Hazard  function 

TriHaz<-function(t){(Tr[TM,l]*tA(Tr[TM,l]- 

l))/(Tr[TM;4]ATr[TM,l])+(Tr[TM,2]*tA(Tr[TM,2]- 

l))/(Tr[TM,5]  ATr[TM,2])+(Tr[TM,3]*tA(Tr[TM,3]-l))/(Tr[TM,6]  ATr[TM,3])} 
####Bi-Weibull  pdf  function 
T  ripdfc-function(t)  {((Tr  [TM,  1]  *t A  (Tr[TM,l]- 
l))/(Tr[TM,4]ATr[TM,l])+(Tr[TM,2]*tA(Tr[TM,2]- 

l))/(Tr[TM,5]ATr[TM,2])+(Tr[TM,3]*tA(Tr[TM,3]-l))/(Tr[TM,6]ATr[TM,3]))*(exp(- 

l*(((t/(Tr[TM,4]))ATr[TM,l])+((t/(Tr[TM,5]))ATr[TM,2])+((t/(Tr[TM,6]))ATr[TM,3 

]))))} 

####Tri-Weibull  CDF 
TWcdf<-function(x,y){ 
bl<-x[l] 
b2<-x[2] 
b3<-x[3] 
tl<-x[4] 
t2<-x[5] 
t3<-x[6] 

F<-c(rep(NA,l)) 

F<-l-(exp(-l*(((y/tl)Abl)+((y/(t2))Ab2)+((y/(t3))Ab3)))) 

return(F) 

} 

TWtestc- 

TWcdf(x=c(Tr[TM,l],Tr[TM,2],Tr[TM,3],Tr[TM,4],Tr[TM,5],Tr[TM,6]),y=Meeker) 


#############  Quad-Weibull  Model  ######################### 

#  qO<-matrix(NA,  nrow=100,  ncol=8,  byrow=FALSE,);  colnames(qO)<- 
c("betal","beta2","beta3","beta4","thetal","theta2","theta3","theta4") 

#  q0[,l]<-c(runif(100,  0.25,5  ))  #betal 

#  q0[,2]<-c(runif(100,  0.25,5  ))  #beta2 

#  q0[,3]<-c(runif(100,  0.25,5  ))  #beta3 

#  q0[,4]<-c(runif(100,  0.25,5  ))  #beta4 

#  q0[,5]<-c(runif(100,  min(Meeker),max(Meeker)))  #thetal 

#  q0[,6]<-c(runif(100,  min(Meeker),max(Meeker)))  #theta2 

#  q0[,7]<-c(runif(100,  min(Meeker),max(Meeker)))  #theta3 
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#  q0[,8]<-c(runif(100,  min(Meeker),max(Meeker)))  #theta4 

# 

#  QQ<-function(x,y,d){ 

#  bl<-x[l] 

#  b2<-x[2] 

#  b3<-x[3] 

#  b4<-x[4] 

#  tl<-x[5] 

#  t2<-x[6] 

#  t3<-x[7] 

#  t4<-x[8] 

# 

#  F<-c(rep(NA,l)) 

# 

#  F<-sumnogffbl*(yAfbl-l))*ftl)A(-bl))+fb2*yAfb2-l))*(ft2)A(-b2))+fb3*yAfb3- 
l))*((t3)A(-b3))+(b4*yACb4-l))*((t4)A(-b4)))*d- 
(C(y)/tl)ACbl)+(y/t2)A[b2)+(y/t3)A(b3)+(y/t4)A(b4))) 

# 

#  return(-F) 

#} 

# 

#  Q<-matrix(NA,  nrow=100,  ncol=9,  byrow=FALSE,);  colnames(Q)<- 
c("betal","beta2","beta3","beta4","thetal","theta2","theta3","theta4", "Value") 

# 

#  for  (i  in  1:100){ 

#  Q[i,]=0 

#  Q[i,l:8]<-optim(par=q0[i,],  fn=QQ,  y=c(Meeker),d=c(Censor),  method="L-BFGS- 
B",  lower=c(0.25,  0.25,  0.25,0.25, 10, 10, 10, 10),  upper=c(115, 115,  90,  90,  500,  500, 
400,  400))  $par 

#  Q[i,9]  <-optim(par=qO[i,],  fn=QQ,  y=c(Meeker),d=c(Censor),  method="L-BFGS- 
B",  lower=c(0.25,  0.25,  0.25,0.25, 10, 10, 10, 10),  upper=c(115, 115,  90,  90,  500,  500, 
400,  400))  $value 

#} 

# 

#  QM<-which.min(Q[,9]) 

#  Q[QM,] 

# 

#  #SEQ<-optim(par=Q[QM,l:8],  fn=QQ,  y=c(Meeker),d=l,  method="L-BFGS-B", 
lower=c(0.25,  0.25,  0.25,0.25, 10, 10, 10, 10),  upper=c(110, 100,  90,  90, 100, 100, 
125, 100),hessian=TRUE)  $hessian 

#  #seq<-sqrt(diag(SEQ)) 

# 

#  ####Quad-Weibull  Reliability  function 
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#  QuadWeib<-function(t){exp(- 

l*(((t/(Q[QM,5]))AQ[QM,l])+((t/(Q[QM,6]))AQ[QM,2])+((t/(Q[QM,7]))AQ[QM,3])+C 

(t/(Q[QM,8]))AQ[QM;4])))} 

#  ####Quad-Weibull  Hazard  function 

#  QuadHaz<-function(t){(Q[QM,l]*tA(Q[QM,l]- 
l))/CQ[QM,5]AQ[QM,l])+(Q[QM,2]*tA(Q[QM,2]- 
l))/(Q[QM,6]AQ[QM,2])+(Q[QM,3]*tA(Q[QM,3]- 

l))/CQ[QM,7]AQ[QM,3])+(Q[QM,4]*tA(Q[QM,4]-l))/(Q[QM;8]AQ[QM,4])} 

#  ####Quad-Weibull  pdf  function 

#  Quadpdf<-function(t){((Q[QM,l]*tA(Q[QM,l]- 
l))/CQ[QM,5]AQ[QM,l])+CQ[QM,2]*tA(Q[QM,2]- 
l))/(Q[QM,6]AQ[QM,2])+(Q[QM,3]*tA(Q[QM,3]- 

l))/(Q[QM,7]AQ[QM,3])+(Q[QM,4]*tA(Q[QM,4]-l))/(Q[QM,8]AQ[QM,4]))*(exp(- 

l*CC(t/CQ[QM,5]))AQ[QM,l])+((t/(Q[QM,6]))AQ[QM,2])+((t/(Q[QM,7]))AQ[QM,3])+C 

(t/(Q[QM,8]))AQ[QM;4]))))} 

# 

#  ####Quad-Weibull  CDF 

#  QWcdf<-function(x,y){ 

#  bl<-x[l] 

#  b2<-x[2] 

#  b3<-x[3] 

#  b4<-x[4] 

#  tl<-x[5] 

#  t2<-x[6] 

#  t3<-x[7] 

#  t4<-x[8] 

# 

#  F<-c(rep(NA,l)) 

# 

#  F<-l-(exp(-l*(((y/tl)Abl)+((y/(t2))Ab2)+((y/(t3))Ab3)+((y/(t4))Ab4)))) 

# 

#  return(F) 

#} 

# 

#  QWtestc- 

QWcdf(x=c(Q[QM,l],Q[QM,2],Q[QM,3],Q[QM,4],Q[QM,5],Q[QM,6],Q[QM,7],Q[QM,8]), 

y=Meeker) 

############  Single  Weibull  Model  ######################### 

Weib<-function(x,y,d){ 

bl<-x[l] 

tl<-x[2] 

n  <  - 1  e  ngth  (Aar  s  e  t2 ) 
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F<-c(rep(NA,l)) 


F<-(d*n*(log(bl)-bl*log(tl))+(d*(bl-l))*sum(log(y))-sum(((y)/tl)Abl)) 

return(-F) 

} 

SW<-matrix(NA,  nrow=100,  ncol=3,  byrow=FALSE,);  colnames(SW)<- 
c("beta", "theta", "Value") 

for  (i  in  1:100){ 

SW[i,]=0 

SW[i,l:2]<-optim(par=c(pO[i,l],  p0[i,3]),  fn=Weib,  y=Meeker,  d=Censor, 
method="L-BFGS-B",  lower=c(0.25,  0.25, 10, 10),  upper=c(5,  5, 1000, 1000))  $par 
SW[i,3]<-Weib(SW[i,l:2],y=Meeker,d=Censor) 

} 

SM<-which.min(SW[,3]) 

####Single-Weibull  Reliability  function 
UniWeib<-function(t){exp(-l*(((t/(SW[SM,2]))ASW[SM,l])))} 
####Single-Weibull  Hazard  function 

UniHaz<-function(t){(SW[SM,l]*tA(SW[SM,l]-l))/(SW[SM,2]ASW[SM,l])} 
####Single-Weibull  pdf  function 

UniPdf<-function(t){(SW[SM,l]*tA(SW[SM,l]-l))/(SW[SM,2]ASW[SM,l])*exp(- 

l*(((t/(SW[SM,2]))ASW[SM,l])))} 

####Single-Weibull  CDF 
SWcdf<-function(x,y){ 
bl<-x[l] 
tl<-x[2] 

F<-c(rep(NA,l)) 

F<-l-exp(-l*(((t/(SW[SM,2]))ASW[SM,l]))) 

return(F) 

} 

SWtest<-SWcdf(x=c(SW[SM,l],SW[SM,2]),y=Meeker) 


#########  Exponentiated  modified  Weibull  extension  Model 
e0<-matrix(NA,  nrow=100,  ncol=4,  byrow=FALSE,);  colnames(p0)<- 
cC'alpha", "lambda", "beta", "gamma") 
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e0[,l]<-c(runif(100,  40,100 ))  #alpha 
e0[,2]<-c(runif(100, 10A-5,10A-2 ))  #lambda 
e0[,3]<-c(runif(100,  .01, .8))  #beta 
e0[,4]<-c(runif(100,  0.25, .625))  #gamma 

EE<-function(x,y,z,d){ 

alp<-x[l] 

lam<-x[2] 

bet<-x[3] 

gam<-x[4] 

n<-length(Meekerl) 

F<-c(rep(NA,l)) 

F<-n*(alp*lam+(l-bet)*log(alp)+log(bet)+log(lam)+log(gam))- 
alp*lam*(sum(exp((y/alp)Abet)))+l/(alpAbet)*(sum(yAbet))+(bet- 
l)*(sum(log(y)))+(gam-l)*sum(l-exp(lam*alp*(l-exp((y/alp)Abet))))+log(prod(l- 
(1-exp  (lam*alp*(l-exp((z/alp)Abet)))Agam))) 

return(-F) 

} 

E<-matrix(NA,  nrow=100,  ncol=5,  byrow=FALSE,);  colnames(P)<- 
c("alpha","lamba", "beta", "gamma",  "Value") 

for  (i  in  1:100){ 

E[i,]=0 

E[i,l:4]<-optim(par=e0[i,],  fn=EE,  y=c(Meekerl),  z=c(Meeker2),d=l,method="L- 
BFGS-B",lower=c(40, 10A-4,  .01,  0.25),  upper=c(100, 10A-2,  .8,  .625))  $par 
E[i,5]<-EE(E[i,l:4],c(Meekerl),c(Meeker2),l) 

> 

EM<-which.min(E[,5]) 

E[EM,] 

####EMWE  Reliability  function 

EWeib<-function(t){l-(l-exp(E[EM,2]*E[EM,l]*(l- 

exp((t/E[EM,l])AE[EM,3]))))A(E[EM,4])} 

####EMWE  Hazard  function 

EHaz<-function(t){(E[EM,2]*E[EM,3]*E[EM,4]*(t/E[EM,l])A(E[EM,3]- 

l)*exp((t/E[EM,l])AE[EM,3]+E[EM,2]*E[EM,l]*(l- 

exp((t/E[EM,l])AE[EM,3]))))/(((l-exp(E[EM,2]*E[EM,l]*(l- 
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exp((t/E[EM,l])AE[EM,3]))))A(l-E[EM,4])+exp(E[EM,2]*E[EM,l]*(l- 

exP((t/E[EM,l])AE[EM,3])))-l))} 

####EMWE  pdf 

EPdf<-function(t){E[EM,2]*E[EM,3]*E[EM,4]*(t/E[EM,l])A(E[EM,3]- 

l)*exp((t/E[EM,l])AE[EM,3]+E[EM,2]*E[EM,l]*(l-exp((t/E[EM,l])AE[EM,3])))*(l- 

exp(E[EM,2]*E[EM,l]*(l-exp((t/E[EM,l])AE[EM,3]))))A(E[EM,4]-l)} 

####EMWE  CDF 

Ecdf<-function(t){(l-exp(E[EM,2]*E[EM,l]*(l- 

exp((t/E[EM,l])AE[EM,3]))))A(E[EM,4])} 

EMtest<-Ecdf(t=eCDF) 

###########################  Plot  Reliability  Functions  against  K-M 
Estimate######################### 

KMFit<-survfit(Surv(Meeker,  Censor)  ~1) 

plot(KMFit,  axes  =  FALSE,  xlab  =  NA,  ylab  =  NA,  xaxs="r",  yaxs="r", 

ylim=range(c(0,l))) 

par(new=TRUE) 

curve(NMWeib,from=0,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab  = 

NA,col="red",lty=3,xaxs="r",  yaxs="r",  ylim=range(c(0,l))) 

par(new=TRUE) 

curve(BiWeib,from=0,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab  = 

NA,xaxs="r",  yaxs="r",  ylim=range(c(0,l))) 

par(new=TRUE) 

curve(TriWeib,from=0,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab  = 

NA,col="orange",xaxs="r",  yaxs="r",  ylim=range(c(0,l))) 

par(new=TRUE) 

curve(EWeib,from=0,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab  = 

NA,col="green",lty=2,xaxs="r",  yaxs="r",  ylim=range(c(0,l))) 

box() 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 
axis(side  =  2,  tck  =  -.015,  labels  =  NA) 
axis(side  =  1,  lwd  =  0,  line  =  -.6) 
axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 
mtext(side  =  1,  "t",  line  =  2) 
mtext(side  =  2,  "R(t)",  line  =  2.5) 

legend("bottomleft",  c("Bi-Weibull","Tri-WeibuH","NMW", 

"EMWE"),lty=c(l,l,3,2),lwd=c(1.5,1.5),col=c("black","orange","red","green"),cex=0.7 

5,bty="n",xjust=0.5) 


###########################  Plot  Density  Functions  against  Censored 
Histogram######################### 

barplot(c(0. 004667, 0.002, 0.002, 0.00 1333, 0.002,0. 002667), space=0,col="white", 
axes=FALSE,  xpd=FALSE,  xaxs="r",yaxs="r",  ylim=range(c(0, 0.005))) 
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par(new=TRUE) 

curve(NMWPdf,from=0,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab  = 
NA,col="red",lty=3,xaxs="r",  yaxs="r",  ylim=range(c(0, 0.005))) 
par(new=TRUE) 

curve(BiPdf,from=0,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab  = 

NA,xaxs="r",  yaxs="r",  ylim=range(c(0, 0.005))) 

par(new=TRUE) 

curve(TriPdf,from=0,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab  = 
NA,col="orange",xaxs="r",  yaxs="r",  ylim=range(c(0, 0.005))) 

#  par(new=TRUE) 

#  curve(EPdf,from=0,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab  = 
NA,col="blue",xaxs="r",  yaxs="r",  ylim=range(c(0, 0.005))) 
par(new=TRUE) 

curve(EPdf,from=0,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab  = 
NA,col="green",lty=2,xaxs="r",  yaxs="r",  ylim=range(c(0, 0.005))) 

#  par(new=TRUE) 

#  curve(UniPdf,from=0,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab  = 
NA,col="Purple",xaxs="r",  yaxs="r",  ylim=range(c(0, 0.005))) 

box() 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1,  "t",  line  =  2) 

mtext(side  =  2,  "f(t)",  line  =  3) 

legend("top",  c("Bi-Weibull","Tri-Weibull","NMW", 

"EMWE"),lty=c(l,l,3,2),lwd=c(1.5,1.5),col=c("black","orange","red","green"),cex=0.7 

5,bty="n",xjust=0.5) 


barplot(c(.005283,  .002791,  .003243,  .0025,  .004444,. 013333),  space=0,col="white", 

axes=FALSE,  xpd=FALSE,  xaxs="i",  ylim=range(c(0,0.015))) 

par(new=TRUE) 

curve(NMWHaz,from=0. 00002,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA, 
ylab  =  NA,col="red",lty=3,xaxs="i",  yaxs="i",  ylim=range(c(0,0.015))) 
par(new=TRUE) 

curve(BiHaz,from=0. 00002,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab 

=  NA,xaxs="i",  yaxs="i",  ylim=range(c(0, 0.015))) 

par(new=TRUE) 

curve(TriHaz,from=0. 00002,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA, 
ylab  =  NA,col="orange",xaxs="i",  yaxs="i",ylim=range(c(0, 0.015))) 

#  par(new=TRUE) 

#  curve(QuadHaz,from=0. 00002,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA, 
ylab  =  NA,col="blue",xaxs="i",  yaxs="i",  ylim=range(c(0, 0.015))) 
par(new=TRUE) 
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curve(EHaz,from=0. 00002,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA,  ylab 
=  NA,col="green",lty=2,xaxs="i",  yaxs="i",  ylim=range(c(0, 0.015))) 

#  par(new=TRUE) 

#  curve(UniHaz,from=0. 00002,  to=max(Meeker),n=1000,  axes  =  FALSE,  xlab  =  NA, 
ylab  =  NA,col="Purple",xaxs="i",  yaxs="i",  ylim=range(c(0,0.015))) 

box() 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1,  "t",  line  =  2) 

mtext(side  =  2,  "h(t)",  line  =  3) 

legend("top",  c("Bi-WeibuH","Tri-Weibull","NMW", 

"EMWE"),lty=c(l,l,3,2),lwd=c(1.5,1.5),col=c("black","orange","red","green"),cex=0.7 

5,bty="n",xjust=0.5) 


N[NM,] 

P[PM,] 

#########  Perform  K-S  Test  for  both  models  and  compare  Against  Meeker  ECDF 

ks.test(RelFreq,  NMtest) 
ks.test(RelFreq,  BWtest) 
ks.test(RelFreq,  TWtest) 
ks.test(RelFreq,  QWtest) 
ks.test(RelFreq,  SWtest) 
ks.test(RelFreq,  EMtest) 


########  Find  standard  errors  for  bi-Weibull  parameter  estimates 

hl<-((P[PM,l]*dataA(P[PM,l]-l))/(P[PM,3]A(P[PM,l]))) 

h2<-((P[PM,2]*dataA(P[PM,2]-l))/(P[PM,4]A(P[PM,2]))) 

rl<-(data/(P[PM,3]))A(P[PM,l]) 

r2<-(data/(P[PM,4]))A(P[PM,2]) 

hh<-((P[PM,l]*dataA(P[PM,l]-l))/(P[PM,3]A(P[PM,l])))+((P[PM,2]*dataA(P[PM,2]- 
1))/(P[PM,4]  A(P[PM,2]))) 

likeblex<-sum(d*(hl/hh*(log(data/P[PM,3])+P[PM,l]A-l)))- 

sum(rl*log(data/P[PM,3])) 

likeb2ex<-sum(d*(h2/hh*(log(data/P[PM,4])+P[PM,2]A-l)))- 

sum(r2*log(data/P[PM,4])) 

liketlex<-sum(-(P[PM,l]/P[PM,3]*d*hl/hh))+sum((P[PM,l]/P[PM,3])*rl) 

Uket2ex<-sum(-(P[PM,2]/P[PM,4]*d*h2/hh))+sum((P[PM,2]/P[PM,4])*r2) 
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likeblex 

likeb2ex 

liketlex 

liket2ex 


Hess<-matrix(NA,nrow=4,  ncol=4,byrow=FALSE) 


Hess[l,l]<-sum((hl*(log(P[PM,3])-log(data))*(log(P[PM,3])-log(data)- 

2/P[PM,l])*hh-(hl*(log(data)-log(P[PM,3]+P[PM,l]A-l)))A2)/(hhA2))- 

sum(rl*(log(data/P[PM,3]))A2) 

Hess[l,2]<-sum((-(hl*(log(data)-log(P[PM,3])+P[PM,l]A-l)*h2*(log(data)- 

log(P[PM,4])+P[PM,2]A-l)))/(hhA2)) 

Hess[l,3]<-sum(((P[PM,l]/P[PM,3])*(-hh*hl*(log(data)-log(P[PM,3])+2*P[PM,l]A 

l)+(P[PM,l]/P[PM,3])*hlA2*(log(data)-log(P[PM,3])+P[PM,l]A- 

l)))/(hhA2))+sum((P[PM,l]/P[PM,3])*rl*(log(data/P[PM,3])+P[PM,l]A-l)) 

Hess[l,4]<-sum((hl*(log(data)-log(P[PM,3])+P[PM,l]A-l)*P[PM,2]*h2*P[PM,4]A- 

l)/(hhA2)) 

Hess[2,l]<-sum((-(h2*(log(data)-log(P[PM,4])+P[PM,2]A-l)*hl*(log(data)- 

log(P[PM,3])+P[PM,l]A-l)))/(hhA2)) 

Hess[2,2]<-sum((h2*(log(P[PM,4])-log(data))*(log(P[PM,4])-log(data)- 

2/P[PM,2])*hh-(h2*(log(data)-log(P[PM,4]+P[PM,2]A-l)))A2)/(hhA2))- 

sum(r2*(log(data/P[PM,4]))A2) 

Hess[2,3]<-sum((h2*(log(data)-log(P[PM,4])+P[PM,2]A-l)*P[PM,l]*hl*P[PM,3]A- 

l)/(hhA2)) 

Hess[2,4]<-sum(((P[PM,2]/P[PM,4])*(-hh*h2*(log(data)-log(P[PM,4])+2*P[PM,2]A 

l)+(P[PM,2]/P[PM,4])*h2A2*(log(data)-log(P[PM,4])+P[PM,2]A- 

l)))/(hhA2))+sum((P[PM,2]/P[PM,4])*r2*(log(data/P[PM,4])+P[PM,2]A-l)) 

Hess[3,l]<-sum(((P[PM,l]/P[PM,3])*(-hh*hl*(log(data)-log(P[PM,3])+2*P[PM,l]A 

l)+(P[PM,l]/P[PM,3])*hlA2*(log(data)-log(P[PM,3])+P[PM,l]A- 

l)))/(hhA2))+sum((P[PM,l]/P[PM,3])*i'l*(log(data/P[PM,3])+P[PM,l]A-l)) 

Hess[3,2]<-sum(((P[PM,l]/P[PM,3])*hl*h2*(log(data)-log(P[PM,4])+P[PM,2]A- 

l))/(hhA2)) 

Hess[3,3]<-sum(((P[PM,l]/P[PM,3])*(P[PM,l]+l)*(P[PM,3]A-l)*hl*hh-(- 

hl*(P[PM,l]/P[PM,3]))A2)/(hhA2))-sum(((P[PM,l]A2+P[PM,l])*rl)/(P[PM,3]A2)) 

Hess[3,4]<-sum(((-P[PM,l]/P[PM,3])*hl*(P[PM,2]/P[PM,4])*h2)/(hhA2)) 

Hess[4,l]<-sum(((P[PM,2]/P[PM,4])*h2*hl*(log(data)-log(P[PM,3])+P[PM,l]A- 

l))/(hhA2)) 

Hess[4,2]<-sum(((P[PM,2]/P[PM,4])*(-hh*h2*(log(data)-log(P[PM,4])+2*P[PM,2]A 
l)+(P[PM,2]/P[PM,4])*h2A2*(log(data)-log(P[PM,4])+P[PM,2]A- 
l)))/(hhA2))+sum((P[PM,2]/P[PM,4])*r2*(log(data/P[PM,4])+P[PM,2]A-l)) 
Hess[4,3]<-sum(((-P[PM,2]/P[PM,4])*h2*(P[PM,l]/P[PM,3])*hl)/(hhA2)) 
Hess[4,4]<-sum(((P[PM,2]/P[PM,4])*(P[PM,2]+l)*(P[PM,4]  A-l)*h2*hh-(- 
h2*(P[PM,2]/P[PM,4]))A2)/(hhA2))-sum(((P[PM,2]A2+P[PM,2])*r2)/(P[PM,4]  A2)) 
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SEc-solve(-Hess) 

SE 

########  Find  standard  errors  for  tri-Weibull  parameter  estimates 
data<-c(PMin  [S  [z]  :N  [z],u]) 

hl<-((Tr[TM,l]*dataA(Tr[TM,l]-l))/(Tr[TM,4]A(Tr[TM,l]))) 

h2<-((Tr[TM,2]*dataA(Tr[TM,2]-l))/(Tr[TM,5]A(Tr[TM,2]))) 

h3<-((Tr[TM,3]*dataA(Tr[TM,3]-l))/(Tr[TM,6]A(Tr[TM,3]))) 

rl<-(data/(Tr[TM,4]))A(Tr[TM,l]) 

r2<-(data/(Tr[TM,5]))A(Tr[TM,2]) 

r3<-(data/(Tr[TM,6]))A(Tr[TM,3]) 

hh<-((Tr[TM,l]*dataA(Tr[TM,l]- 

l))/(Tr[TM,4]A(Tr[TM,l])))+((Tr[TM,2]*dataA(Tr[TM,2]- 
l))/(Tr[TM,5]A(Tr[TM,2])))+((Tr[TM,3]*dataA(Tr[TM,3]- 
1)) /(Tr[TM,6]  A(Tr[TM,3]))) 

likeblex<-smn((hl*(log(data)-log(Tr[TM,4])+Tr[TM,l]A-l))/hh)- 

sum(rl*log(data/Tr[TM,4])) 

likeb2ex<-sum((h2*(log(data)-log(Tr[TM,5])+Tr[TM,2]A-l))/hh)- 

sum(r2*log(data/Tr[TM,5])) 

likeb3ex<-sum((h3*(log(data)-log(Tr[TM,6])+Tr[TM,3]A-l))/hh)- 

sum(r3*log(data/Tr[TM,6])) 

liketlex<-sum(-(Tr[TM,l]/Tr[TM,4]*hl/hh))+sum((Tr[TM,l]/Tr[TM,4])*rl) 

liket2ex<-sum(-(Tr[TM,2]/Tr[TM,5]*h2/hh))+sum((Tr[TM,2]/Tr[TM,5])*r2) 

liket3ex<-sum(-(Tr[TM,3]/Tr[TM,6]*h3/hh))+sum((Tr[TM,3]/Tr[TM,6])*r3) 

likeblex 

likeb2ex 

likeb3ex 

liketlex 

liket2ex 

liket3ex 

Hess<-matrix(NA,nrow=6,  ncol=6,byrow=FALSE) 

Hess[l,l]<-sum((hl*(log(Tr[TM,4])-log(data))*(log(Tr[TM,4])-log(data)- 

2/Tr[TM,l])*hh-(hl*(log(data)-log(Tr[TM,4]+Tr[TM,l]A-l)))A2)/(hhA2))- 

sum(rl*(log(data/Tr[TM,4]))A2) 

Hess[l,2]<-sum((-(hl*(log(data)-log(Tr[TM,4])+Tr[TM,l]A-l)*h2*(log(data) 

log(Tr[TM,5])+Tr[TM,2]A-l)))/(hhA2)) 

Hess[l,3]<-sum((-(hl*(log(data)-log(Tr[TM;4])+Tr[TM,l]A-l)*h3*(log(data) 

log(Tr[TM,6])+Tr[TM,3]A-l)))/(hhA2)) 
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Hess[l,4]<-sum(((Tr[TM,l]/Tr[TM,4])*(-hh*hl*(log(data)- 

log(Tr[TM;4])+2*Tr[TM,l]A-l)+(Tr[TM,l]/Tr[TM;4])*hlA2*(log(data)- 

log(Tr[TM,4])+Tr[TM,l]A- 

l)))/(hhA2))+sum((Tr[TM,l]/Tr[TM,4])*rl*(log(data/Tr[TM,4])+Tr[TM,l]A-l)) 
Hess[l,5]<-sum((hl*(log(data)4og(Tr[TM,4])+Tr[TM,l]A- 
l)*Tr[TM,2]*h2*Tr[TM,5]A-l)/(hhA2)) 
Hess[l,6]<-sum((hl*(log(data)4og(Tr[TM,4])+Tr[TM,l]A- 
l)*Tr[TM,3]*h3*Tr[TM,6]  A-l)/(hhA2)) 

Hess[2,l]<-sum((-(h2*(log(data)4og(Tr[TM,5])+Tr[TM,2]A-l)*hl*(log(data)- 

log(Tr[TM;4])+Tr[TM,l]A-l)))/(hhA2)) 

Hess[2,2]<-sum((h2*(log(Tr[TM,5])-log(data))*(log(Tr[TM,5])-log(data)- 

2/Tr[TM,2])*hh-(h2*Qog(data)4og(Tr[TM,5]+Tr[TM,2]A-l)))A2)/(hhA2))- 

sum(r2*(log(data/Tr[TM,5]))A2) 

Hess[2,3]<-sum((-(h2*(log(data)-log(Tr[TM,5])+Tr[TM,2]A-l)*h3*(log(data)- 

log(Tr[TM;6])+Tr[TM,3]A-l)))/(hhA2)) 

Hess[2,4]<-sum((h2*(log(data)-log(Tr[TM,5])+Tr[TM,2]A- 

l)*Tr[TM,l]*hl*Tr[TM,4]A-l)/(hhA2)) 

Hess[2,5]<-sum(((Tr[TM,2]/Tr[TM,5])*(-hh*h2*(log(data)- 

log(Tr[TM,5])+2*Tr[TM,2]A-l)+(Tr[TM,2]/Tr[TM,5])*h2A2*(log(data)- 

log(Tr[TM,5])+Tr[TM,2]A- 

l)))/(hhA2))+sum((Tr[TM,2]/Tr[TM,5])*r2*(log(data/Tr[TM,5])+Tr[TM,2]A-l)) 
Hess[2,6]<-sum((h2*(log(data)-log(Tr[TM,5])+Tr[TM,2]A- 
l)*Tr[TM,3]*h3*Tr[TM,6]  A-l)/(hhA2)) 

Hess[3,l]<-sum((-(h3*(log(data)-log(Tr[TM,6])+Tr[TM,3]A-l)*hl*(log(data)- 

log(Tr[TM,4])+Tr[TM,l]A-l)))/(hhA2)) 

Hess[3,2]<-sum((-(h3*(log(data)-log(Tr[TM,6])+Tr[TM,3]A-l)*h2*(log(data)- 

log(Tr[TM,5])+Tr[TM,2]A-l)))/(hhA2)) 

Hess[3,3]<-sum((h3*(log(Tr[TM,6])-log(data))*(log(Tr[TM;6])-log(data)- 

2/Tr[TM,3])*hh-(h3*(log(data)-log(Tr[TM,6]+Tr[TM,3]A-l)))A2)/(hhA2))- 

sum(r3*(log(data/Tr[TM,6]))A2) 

Hess[3,4]<-sum((h3*(log(data)-log(Tr[TM,6])+Tr[TM,3]A- 

l)*Tr[TM,l]*hl*Tr[TM,4]A-l)/(hhA2)) 

Hess[3,5]<-sum((h3*(log(data)4og(Tr[TM,6])+Tr[TM,3]A- 

l)*Tr[TM,2]*h2*Tr[TM,5]A-l)/(hhA2)) 

Hess[3,6]<-sum(((Tr[TM,3]/Tr[TM,6])*(-hh*h3*(log(data)- 

log(Tr[TM,6])+2*Tr[TM,3]A-l)+(Tr[TM,3]/Tr[TM,6])*h3A2*(log(data)- 

log(Tr[TM,6])+Tr[TM,3]A- 

l)))/(hhA2))+sum((Tr[TM,3]/Tr[TM,6])*r3*(log(data/Tr[TM;6])+Tr[TM,3]A-l)) 

Hess[4;l]<-sum(((Tr[TM,l]/Tr[TM;4])*(-hh*hl*(log(data)- 

log(Tr[TM,4])+2*Tr[TM,l]A-l)+(Tr[TM,l]/Tr[TM;4])*hlA2*(log(data)- 

log(Tr[TM,4])+Tr[TM,l]A- 

l)))/(hhA2))+sum((Tr[TM,l]/Tr[TM,4])*rl*(log(data/Tr[TM,4])+Tr[TM,l]A-l)) 

Hess[4,2]<-sum(((Tr[TM4]/Tr[TM,4])*hl*h2*(log(data)-log(Tr[TM,5])+Tr[TM,2]A- 

l))/(hhA2)) 
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Hess[4,3]<-sum(((Tr[TM,l]/Tr[TM,4])*hl*h3*(log(data)-log(Tr[TM,6])+Tr[TM,3]A 

l))/(hhA2)) 

Hess[4,4]<-sum(((Tr[TM,l]/Tr[TM,4])*(Tr[TM,l]+l)*(Tr[TM,4]A-l)*hl*hh-(- 

hl*(Tr[TM,l]/Tr[TM,4]))A2)/(hhA2))- 

sum(((Tr[TM,l]A2+Tr[TM,l])*rl)/(Tr[TM,4]A2)) 

Hess[4;5]<-sum(((-Tr[TM,l]/Tr[TM;4])*hl*(Tr[TM,2]/Tr[TM,5])*h2)/(hhA2)) 

Hess[4,6]<-sum(((-Tr[TM,l]/Tr[TM;4])*hl*(Tr[TM,3]/Tr[TM,6])*h3)/(hhA2)) 

Hess[5,l]<-sum(((Tr[TM,2]/Tr[TM,5])*h2*hl*(log(data)-log(Tr[TM,4])+Tr[TM,l]A 

l))/(hhA2)) 

Hess[5,2]<-sum(((Tr[TM,2]/Tr[TM,5])*(-hh*h2*(log(data)- 

log(Tr[TM,5])+2*Tr[TM,2]A-l)+(Tr[TM,2]/Tr[TM,5])*h2A2*(log(data)- 

log(Tr[TM,5])+Tr[TM,2]A- 

l)))/(hhA2))+sum((Tr[TM,2]/Tr[TM,5])*r2*(log(data/Tr[TM,5])+Tr[TM,2]A-l)) 

Hess[53]<-sum(((Tr[TM,2]/Tr[TM,5])*h2*h3*(log(data)-log(Tr[TM,6])+Tr[TM,3]A 

l))/(hhA2)) 

Hess[5,4]<-sum(((-Tr[TM,2]/Tr[TM,5])*h2*(Tr[TM,l]/Tr[TM,4])*hl)/(hhA2)) 
Hess[5,5]<-sum(((Tr[TM,2]/Tr[TM,5])*(Tr[TM,2]+l)*(Tr[TM,5]A-l)*h2*hh-(- 
h2*(Tr[TM,2]/Tr[TM,5]))A2)/(hhA2))- 
sum(((Tr[TM,2]A2+Tr[TM,2])*r2) /(Tr  [TM,5]  A2)) 

Hess[5,6]<-sum(((-Tr[TM,2]/Tr[TM,5])*h2*(Tr[TM,3]/Tr[TM,6])*h3)/(hhA2)) 

Hess[6,l]<-sum(((Tr[TM3]/Tr[TM,6])*h3*hl*(log(data)-log(Tr[TM,4])+Tr[TM,l]A 

l))/(hhA2)) 

Hess[6,2]<-sum(((Tr[TM,3]/Tr[TM,6])*h3*h2*(log(data)-log(Tr[TM,5])+Tr[TM,2]A 

l))/(hhA2))  . 

Hess[6,3]<-sum(((Tr[TM,3]/Tr[TM;6])*(-hh*h3*(log(data)- 

log(Tr[TM,6])+2*Tr[TM,3]A-l)+(Tr[TM,3]/Tr[TM,6])*h3A2*(log(data)- 

log(Tr[TM,6])+Tr[TM,3]A- 

l)))/(hhA2))+sum((Tr[TM,3]/Tr[TM,6])*r3*(log(data/Tr[TM;6])+Tr[TM,3]A-l)) 
Hess[6,4]<-sum(((-Tr[TM3]/Tr[TM,6])*h3*(Tr[TM,l]/Tr[TM,4])*hl)/(hhA2)) 
Hess[6,5]<-sum(((-Tr[TM3]/Tr[TM,6])*h3*(Tr[TM,2]/Tr[TM,5])*h2)/(hhA2)) 
Hess[6;6]<-sum(((Tr[TM,3]/Tr[TM;6])*(Tr[TM,3]  +  l)*(Tr[TM,6]A-l)*h3*hh-(- 
h3*(Tr[TM,3]/Tr[TM,6]))A2)/(hhA2))- 
sum(((Tr[TM,3]A2+Tr[TM,3])*r3)/(Tr[TM,6]  A2)) 


SE<-solve(-Hess) 

SE 
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Appendix  D:  Computer  Code  Used  to  Simulate  and  Analyze  Poly-Weibull 
Failure  Observations  as  Discussed  in  Chapter  V 


###  set.seed(42) 

###  This  code  generates  observations  from  n  independent,  distinct  Weibull 
###  distributions  denoting  the  n  independent  failure  modes  in  each  system.  The 
###  minimum  of  these  observations  is  kept  as  the  actual  system  failure  time,  n  is 
###  simulated  as  either  10,20,  40  or  100.  The  failure  data  thus  comes  from  a  Poly- 
###  Weibull  distribution.  The  code  captures  this  data  into  a  matrix  PMin  for 
###  various  combinations  of  number  of  systems  and  number  of  failure  modes 

########  CONTROL  VARIABLES 

a<-5  ##  Selects  the  ath  element  of  vector  z  for  the  number  of  systems 

starts<-10  ##  The  number  of  starting  points  in  the  parameter  space 
quality<-10  ##  Value  for  the  shape  parameter  theta 
f<-10000  ##  The  number  of  iterations 

z<-as.vector(c(10,20,40,100))  #Observations 

########  Initialize  matrices  of  parameters  values  for  each  model 
SW<-matrix(NA,  nrow=starts,  ncol=3,  byrow=FALSE,) 

P<-matrix(NA,  nrow=starts,  ncol=5,  byrow=FALSE,) 

Tr<-matrix(NA,  nrow=starts,  ncol=7,  byrow=FALSE,) 

PROBs  <-matrix(NA,  nrow=f,  ncol=5,  byrow=TRUE) 

PROBb  <-matrix(NA,  nrow=f,  ncol=5,  byrow=TRUE) 

PROBt  <-matrix(NA,  nrow=f,  ncol=5,  byrow=TRUE) 

OSLM  <-matrix(NA,  nrow=f,  ncol=5,  byrow=TRUE) 

OSL5M  <-matrix(NA,  nrow=f,  ncol=5,  byrow=TRUE) 

OSL10M<-matrix(NA,  nrow=f,  ncol=5,  byrow=TRUE) 

OSL20M<-matrix(NA,  nrow=f,  ncol=5,  byrow=TRUE) 

MEAN<-matrix(NA,  nrow=7,  ncol=5,  byrow=TRUE) 

colnames(P)  <-c("betal","beta2","thetal","theta2",  "Value") 
colnames(SW)<-c("beta", "theta", "Value") 

colnames(Tr)<-c("betal","beta2","beta3","thetal","theta2","theta3", "Value") 
colnames(PROBb)  <-c("5-modes","  10-modes", "20-modes", "40-modes", "50-modes") 
colnames(PROBs)  <-c("5-modes","  10-modes", "20-modes", "40-modes", "50-modes") 
colnames(PROBt)  <-c("5-modes","  10-modes", "20-modes", "40-modes", "50-modes") 
colnames(OSLM)<-c("5-modes", "10-modes", "20-modes", "40-modes", "50-modes") 
colnames(OSL5M)<-c("5-modes","  10-modes", "20-modes", "40-modes", "50-modes") 
colnames(OSL10M)<-c("5-modes","  10-modes", "20-modes", "40-modes", "50- 
modes") 
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colnames(OSL20M)<-c("5-modes","  10-modes",  "20-modes", "40-modes", "50- 
modes") 

rownames(MEAN)<-c("Weibull",  "bi-Weibull",  "tri-Weibull","AD(p- 
value)","OSL5","OSL10","OSL20") 

colnames(MEAN)<-c("5  modes'V'10  modes", "20  modes", "40  modes", "50  modes") 
for  (c  in  l:f)  { 

AIC<-matrix(NA,  nrow=5,  ncol=ll,  byrow=TRUE);  colnames(AIC)<-c("Weib","Bi- 
Weib","Tri-Weib","AIC","p-S","p-B","p-T","AD(p-value)",  "OSL5","OSL10","OSL20") 

###########################  Single  Weibull  Model 
######################### 

sO<-matrix(NA,  nrow=starts,  ncol=2,  byrow=FALSE,);  colnames(sO)<-c("a","b") 
sO[,l]<-c(runif(starts,  .25,10))  #beta  -  shape  parameter 
s0[,2]<-c(runif(starts,  1,1500))  #alpha  -  scale  parameter 

###########################  Bi-Weibull  Model 
######################### 

p0<-matrix(NA,  nrow=starts,  ncol=4,  byrow=FALSE,);  colnames(p0)<- 

c("betal","beta2","thetal","theta2") 

p0[,l]<-c(runif(starts,  .15, .9  ))  #betal 

p0[,2]<-c(runif(starts,  5.0,40 ))  #beta2 

p0[,3]<-c(runif(starts,  2.5,1000))  #thetal 

p0[,4]<-c(runif(starts,  2.5,1000))  #theta2 

########  Tri-Weibull  Model 

t0<-matrix(NA,  nrow=starts,  ncol=6,  byrow=FALSE,) 

colnames(t0)<-c("betal","beta2","beta3","thetal","theta2","theta3") 

t0[,l]<-c(runif(starts,  .15, .9))  #betal 
t0[,2]<-c(runif(starts,  5,30 ))  #beta2 
t0[,3]<-c(runif(starts,  1.25,40 ))  #beta3 
t0[,4]<-c(runif(starts,  2.5,1000))  #thetal 
t0[,5]<-c(runif(starts,  2.5,1000))  #theta2 
t0[,6]<-c(runif(starts,  2.5,1000))  #theta3 


j<-l 

i<-l 


######  INITIALIZE  THE  MATRIX  THAT  THE  FAILURE  OBSERVATIONS  WILL  BE 
WRITTEN  TO 

PMin<-matrix(NA,  nrow=z[a],  ncol=5,  byrow=FALSE,) 

colnames(PMin)<-c("5-Modes","  10-Modes", "20-Modes",1 "40-Modes",1 "50-Modes") 
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m<-as.vector(c(5, 10,  20,  40,  50)) 

for  (j  in  1:5)  {  #  grab  an  element  from  the  vector  m 

#######  Create  Weibull  shape  and  scale  parameters  for  each  failure  mode 
CV<-matrix(NA,  nrow=m[j],  ncol=2,  byrow=FALSE);  colnames(CV)<-c("CoV", 
"Mean") 

NLM<-matrix(NA,  nrow=m[j],  ncol=2,  byrow=FALSE);  colnames(NLM)<-c("Betaj", 
"Alphaj") 

BetaP<-as.matrix(c(1.25,quality)) 

#########  Simulate  a  coefficient  of  variation  (COV)  value 
CV[,l]<-c(rbeta(m[j],  BetaP[l,],  BetaP[2,])*(5.5-.05)+.05) 

#########  Define  Mean 

CV[,2]<-c(rbeta(m[j],  (l+((5.5-CV[,l])/(5.5-.05)))A2,  (l+((CV[,l]-.05)/(5.5- 
.05)))A2)*3000) 

##############Functions############### 

COV<-function(x,y)  |betaj<-x[l] 

(((((gamma(l+2/betaj))/(gamma(l+l/betaj))A2)-l)A0.5)-y) 

} 

#########  Generate  Solutions  ############# 
for  (i  in  l:m[j])  { 

NLM[i,l]<-uniroot(COV,  c(0.1,  65),  y=c(CV[i,l]))  $root 
NLM  [i,2]  <-CV[i,2]  /  (gamma(l+ 1/NLM  [i,  1])) 

} 

b<-c(NLM[,l]) 

t<-c(NLM[,2]) 

for  (i  in  l:z[a])  { 

###Generate  possible  failure  time  for  each  failure  mode  and  select  the  minimum 
value 

###as  the  observed  system  failure  time 
PMin[i,j]<-min(rweibull(m[j],  b,  t)) 


} 

data<-sort(PMin[,j]) 
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#########  Compare  Fit  of  the  Weibull,  Poly-Weibull  and  New  Modified  Weibull 
Models 

###########################  Single  Weibull  Model 
######################### 

Weib<-function(x,y){bl<-x[l];  tl<-x[2] 

F<-c(rep(NA,l)) 

F<-(z[a]*(log(bl)-bl*log(tl))+(bl-l)*sum(log(y))-sum(((y)/tl)Abl)) 

return(-F) 

} 

SW<-matrix(NA,  nrow=starts,  ncol=3,  byrow=FALSE,);  colnames(SW)<- 
c("beta", "theta", "Value") 

for  (i  in  1: starts)! 

SW[i,]=0 

SW[i,l:2]<-nlminb(c(sO[i,]),  Weib,  y=c(PMin[,j]),  lower=c(.25,.5), 
upper=c(40,750))  $par 

SW[i,3]<-Weib(SW[i,l:2],c(PMin[,j]))  #value  of  the  negative  log  likelihood 
function 
} 

SM<-which.min(SW[,3])###  Select  the  minimum  value  of  the  10  solutions 

###Perform  the  Anderson-Darling  Goodness  of  Fit  test  for  each  data  set 
AD<-matrix(NA,  nrow=z[a],  ncol=l,  byrow=TRUE) 
for  (i  in  l:z[a])  { 

###  Weibull  Modified  test  statistic 

AD  [i,l]<-(log(l-exp(-l*(data[i]/SW[SM,2])A(SW[SM,l])))-((data[z[a]+l- 
i])/SW[SM,2])A(SW[SM,l]))*((l-2*i)/z[a]) 

} 

ADstar<-(sum(AD  [,l])-z[a])*(l+0.2/sqrt(z[a])) 

()SL<-0 

if(ADstar<.474)  OSL<-.25 

if(.474<=ADstar&&ADstar<=  .637)  OSL<-(.l-.25)/(.637-.474)*(ADstar-.474)+.25 
if(.637<=ADstar&&ADstar<=  .757)  OSL<-(.05-.l)/(.757-.637)*(ADstar-.637)+.10 
if(.757<=ADstar&&ADstar<=  .877)  OSL<-(.025-.05)/(.877-.757)*(ADstar-.757)+.05 
if(.877<=ADstar&&ADstar<=1.038)  OSL<-(.01-.025)/(1.038-.877)*(ADstar- 
.877)+. 025 

if(ADstar>1.038)  OSL<-.01 

#######  Perform  Bernoulli  test  for  each  data  set 
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0SL5  <-ifelse(OSL>0.05, 1,  0) 

OSL10<-ifelse(OSL>0.10, 1,  0) 

OSL20<-ifelse(OSL>0.20, 1,  0) 

#######  Compute  the  corrected  AIC  value  for  the  single  Weibull 
AIC[j,l]<-2*2-2*(-SW[SM,3])+((2*2*(2+l))/(z[a]-2-l)) 

#######  Bi-Weibull  Model 
GG<-function(x,y,d){ 
bl<-x[l];  b2<-x[2];  tl<-x[3];  t2<-x[4] 

F<-c(rep(NA,l)) 

F<-sum(log(((bl*(yA(bl-l))*(tl)A(-bl))+(b2*yA(b2-l))*((t2)A(-b2))))*d- 

(((y)/tl)ACbl)+(y/t2)A[b2))) 

return(-F) 

} 

P<-matrix(NA,  nrow=starts,  ncol=5,  byrow=FALSE,);  colnames(P)<- 
c("betal","beta2","thetal","theta2",  "Value") 

for  (i  in  1: starts)! 

P[iJ=0 

P[i,l:4]<-nlminb(c(p0[i,l:4]),  GG,  y=c(PMin[,j]),d=l,  lower=c(.25,.25,.5,.5), 
upper=c(70, 100, 1000, 1000))  $par 
P[i,5]<-GG(P[i,l:4],c(PMin[,j]),l) 

} 

PM<-which.min(P[,5]);P[PM,] 

AIC[j,2]<-2*4-2*(-P[PM,5])  +((2*4*(4+l))/(z[a]-4-l)) 

###########################  Tri-Weibull  Model 
######################### 

TT  <-function(x,y,d){ 

bl<-x[l];  b2<-x[2];  b3<-x[3];  tl<-x[4];  t2<-x[5];  t3<-x[6] 

F<-c(rep(NA,l)) 

F<-sum(logffbl*(yAfbl-l))*(tl)Af-bl))  +  fb2*yAfb2-l))*((t2)A(-b2))+fb3*yAfb3- 
l))*C(t3)A(-b3)))*d-(((y)/tl)A(bl)+(y/t2)A(b2)+(y/t3)A(b3))) 

return(-F) 

} 
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Tr<-matrix(NA,  nrow=starts,  ncol=7,  byrow=FALSE,);  colnames(Tr)<- 
c("betal","beta2","beta3","thetal","theta2","theta3", "Value") 

for  (i  in  1: starts)! 

Tr[i,]=0 

Tr[i,l:6]<-nlminb(c(t0[i,l:6]),  TT,y=c(PMin[,j]),d=l,  lower=c(.15,.15,.15,.5,.5,.5), 
upper=c(100, 100, 100, 1000, 1000, 1000))  $par 

Tr[i,7]<-TT(Tr[U:6],c(PMintj]),l) 

} 

TM<-which.min(Tr[,7]);Tr[TM,] 

AIC[j,3]  <-2*6-2*(-Tr[TM,7])+((2*6*(6+l))/(z[a]-6-l)) 

AIC[j,4]  <-which.min(AIC[j,l:3]) 

AIC[j,5]  <-exp((AIC[j,AIC[j,4]]-AIC[j,l])/2) 

AIC [j,6]  <-exp((AIC[j,AIC[j,4]]-AIC[j,2])/2) 

AIC [j,7]  <-exp((AIC[j,AIC[j,4]]-AIC[j,3])/2) 

AIC[j,8]  <-OSL 
AIC [j,9]  <-0SL5 
AIC [j,10]  <-OSL10 
AIC [j, 11]  <-OSL20 

} 


PROBs  [c,]  <-AIC[,5] 

PROBb  [c,]  <-AIC[,6  ] 

PROBt  [c,]  <-AIC[,7  ] 

OSLM  [c,]  <-AIC[,8] 

OSL5M  [c,]  <-AIC[,9  ] 

OSL10M  [c,]  <-AIC[,10] 

OSL20M  [cj  <-AIC[,ll] 

} 

MEAN[l,]<-c(mean(PROBs  [,l]),mean(PROBs  [,2]),mean(PROBs  [,3]),mean(PROBs 

[,4]),mean(PROBs  [,5])) 

MEAN[2,]<-c(mean(PROBb  [,l]),mean(PROBb  [,2]),mean(PROBb 

[,3]),mean(PROBb  [,4]),mean(PR0Bb  [,5])) 

MEAN[3,]<-c(mean(PROBt  [,l]),mean(PROBt  [,2]),mean(PROBt  [,3]),mean(PROBt 

[,4]),mean(PR0Bt  [,5])) 

MEAN[4,]<-c(mean(0SLM  [,l]),mean(OSLM  [,2]),mean(OSLM  [,3]),mean(OSLM 

[,4]),mean(0SLM  [,5])) 

MEAN[5,]<-c(sum(OSL5M  [,l])/f,sum(OSL5M  [,2])/f,sum(OSL5M 

[,3])/f,sum(OSL5M  [,4])/f,sum(OSL5M  [,5])/f) 

MEAN[6,]<-c(sum(OSL10M 

[,l])/f,sum(OSL10M[,2])/f,sum(OSL10M[,3])/f,sum(OSL10M[,4])/f,sum(OSL10M[,5 

m 
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MEAN[7,]<-c(sum(OSL20M 

[,l])/f,sum(OSL20M[,2])/f,sum(OSL20M[,3])/f,sum(OSL20M[,4])/f,sum(OSL20M[,5 

m 

###Diagnostic  plots  to  verify  assess  the  code 

qqPlot(c(PMin[,j]),  dist="weibull",  shape=SW[SM,l],  scale=SW[SM,2]) 
CCCC<-sort(c(PMin[,j])) 

Ppos<-pp(CCCC,  a=0) 

quant<-(log(l/(l-Ppos))) 

plot(loglO(c(CCCC))~loglO(quant)) 

UniPdf<-function(t){(SW[SM,l]*tA(SW[SM,l]-l))/(SW[SM,2]ASW[SM,l])*exp(- 

l*(((t/(SW[SM,2]))ASW[SM,l])))} 

BiPdf<-function(t){((P[PM,l]*tA(P[PM,l]- 

l))/(P[PM,3]AP[PM,l])+(P[PM,2]*tA(P[PM,2]-l))/(P[PM;4]AP[PM,2]))*(exp(- 

l*(C(t/(P[PM,3]))AP[PM,l])+((t/(P[PM;4]))AP[PM,2]))))} 

TriPdf<-function(t){((Tr[TM,l]*tA(Tr[TM,l]- 

l))/(Tr[TM;4]ATr[TM,l])+(Tr[TM,2]*tA(Tr[TM,2]- 

l))/(Tr[TM,5]  ATr[TM,2])+(Tr[TM,3]*tA(Tr[TM,3]-l))/(Tr[TM,6]  ATr[TM,3]))*(exp(- 
l*(((V(Tr[TM,4]))ATr[TM,l])+((t/(Tr[TM,5]))ATr[TM,2])+((t/(Tr[TM,6]))ATr[TM,3 
]))))} 

hist(c(PMin[,j]))#,  breaks=z[a]/4) 
par(new=TRUE) 

curve(UniPdf,  from=2.5-min(PMin[,j]),to=max(PMin[,j]),  col="red",n=2000,  xlab="", 

ylab="",  axes=FALSE,  xaxs="r",  yaxs="r") 

par(new=TRUE) 

curve(BiPdf,  from=2.5-min(PMin[,j]),to=max(PMin[,j]),  col="blue",n=2000,  xlab="", 

ylab="",  axes=FALSE,  xaxs="r",  yaxs="r") 

par(new=TRUE) 

curve(TriPdf,  from=2.5-min(PMin[,j]),to=max(PMin[,j]),  col="green",n=2000, 
xlab="",  ylab="",  axes=FALSE,  xaxs="r",  yaxs="r") 

KMFit<-survfit(Surv(c(PMin[,j]),  c(rep(l,z[a])))~l) 

#####Single-Weibull  Reliability  function 

UniWeib<-function(time){exp(-l*(((((time))/(SW[SM,2]))ASW[SM,l])))} 

#####Bi-Weibull  Reliability  function 
BiWeib<-function(time) 

{exp(-l*((((time)/(P[PM,3]))AP[PM,l])+(((time)/(P[PM,4]))AP[PM,2])))} 
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#####Tri-Weibull  Reliability  function 
TriWeib<-function(time) 

{exp(-l*((((time)/(Tr[TM,4]))ATr[TM,l])+(((time)/(Tr[TM,5]))ATr[TM,2])+ 

C((time)/(Tr[TM,6]))ATr[TM,3])))} 

########  Plot  Reliability  Functions  against  K-M  Estimate 
plot(KMFit,  axes  =  FALSE,  xlab  =  NA,  ylab  =  NA,  xaxs="r",  yaxs="r", 
ylim=range(c(0,l))) 
par(new=TRUE) 

curve(UniWeib,from=0,to=max(PMin[,j]),n=100,axes=FALSE,xlab=NA,ylab=NA,col=" 

red",lty=3,xaxs="r",yaxs="r",ylim=range(c(0,l))) 

par(new=TRUE) 

curve(BiWeib,from=0,to=max(PMin[,j]),n=100,axes=FALSE,xlab=NA,ylab=NA,lty=l, 

xaxs="r",yaxs="r",  yiim=range(c(0,l))) 

par(new=TRUE) 

curve(TriWeib,from=0,  to=max(PMin[,j]),n=100,  axes  =  FALSE,  xlab  =  NA, 
ylab=NA,col="orange",xaxs="r",  yaxs="r",  yiim=range(c(0,l))) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 
axis(side  =  2,  tck  =  -.015,  labels  =  NA) 
axis(side  =  1,  lwd  =  0,  line  =  -.6) 
axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 
mtext(side  =  l,"t",  line  =  2.0) 
mtext(side  =  2,  "R(t)",  line  =  2.5) 


##  Generate  Histograms 

g<-c(0, . 05, .1,. 15, .2, .25, .3, .35, .4, .45, .5, .55, .6, .65, .7, .75, .8, .85, .9, .95,1) 
setwd("C:/Users/Jason/Desktop") 


# 

# 


tiff(filename  =  "5byl00.tiff", 

width  =  1800,  height  =  2700,  units  =  "px",  pointsize  =  12, 
compression  =  c("none"),bg  =  "white",  res  =  300) 


par(mfrow  =  c(3,l),oma=rep(4,  4),  mar=rep(l,  4) ) 


hist(PROBs[,l],  col=rgb(l,0,0,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),  ylim=c(0,f),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 
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mtext(side  =  2,  "frequency",  line  =  2.5) 
mtext(side  =  4,  "Weibull",  line  =  .5) 

hist(PROBb[,l],  col=rgb(0,l,0,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),  ylim=c(0,f),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 

mtext(side  =  4,  "Bi-Weibull",  line  =  .5) 

hist(PROBt[,l],col=rgb(0,0,l,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),ylim=c(0,f),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 

mtext(side  =  4,  "Tri-Weibull",  line  =  .5) 

mtext(  expression("AIC"[corrected]*"  (5, 100, Low)" ),  outer  =  TRUE) 
dev.off(2) 

# 

# 

tiff(filename  =  "lObylOO.tiff", 

width  =  1800,  height  =  2700,  units  =  "px",  pointsize  =  12, 
compression  =  c("none"),bg  =  "white",  res  =  300) 

par(mfrow  =  c(3,l),oma=rep(4,  4),  mar=rep(l,  4) ) 

hist(PROBs[,2],  col=rgb(l,0,0,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 

mtext(side  =  4,  "Weibull",  line  =  .5) 
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hist(PR0Bb[,2],  col=rgb(0,l,0,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 

mtext(side  =  4,  "Bi-Weibull",  line  =  .5) 

hist(PROBt[,2],  col=rgb(0,0,l,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 

mtext(side  =  4,  "Tri-Weibull",  line  =  .5) 

mtext(  expression("AIC"[corrected]*"  (10, 100, Low)" ),  outer  =  TRUE) 
dev.off(2) 

# 

# 

tiff(filename  =  "20byl00.tiff", 

width  =  1800,  height  =  2700,  units  =  "px",  pointsize  =  12, 
compression  =  c("none"),bg  =  "white",  res  =  300) 

par(mfrow  =  c(3,l),oma=rep(4,  4),  mar=rep(l,  4) ) 

hist(PROBs[,3],col=rgb(l,0,0,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, "",  line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 

mtext(side  =  4,  "Weibull",  line  =  .5) 

hist(PROBb[,3],  col=rgb(0,l,0,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 
xaxs="r",  yaxs="r",  xlim=c(0,l),  breaks=g) 
axis(side  =  1,  tck  =  -.015,  labels  =  NA) 
axis(side  =  2,  tck  =  -.015,  labels  =  NA) 
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axis(side  =  1,  lwd  =  0,  line  =  -.6) 
axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 
mtext(side  =  1, line  =  2.0) 
mtext(side  =  2,  "frequency",  line  =  2.5) 
mtext(side  =  4,  "Bi-Weibull",  line  =  .5) 

hist(PROBt[,3],col=rgb(0,0,l,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 

mtext(side  =  4,  "Tri-Weibull",  line  =  .5) 

mtext(  expression("AIC"[corrected]*"  (20, 100, Low)" ),  outer  =  TRUE) 
dev.off(2) 

# 

# 

tiff(filename  =  "40byl00.tiff", 

width  =  1800,  height  =  2700,  units  =  "px",  pointsize  =  12, 
compression  =  c("none"),bg  =  "white",  res  =  300) 

par(mfrow  =  c(3,l),oma=rep(4,  4),  mar=rep(l,  4) ) 

hist(PROBs[,4],col=rgb(l,0,0,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 

mtext(side  =  4,  "Weibull",  line  =  .5) 

hist(PR0Bb[,4],  col=rgb(0,l,0,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 
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mtext(side  =  4,  "Bi-Weibull",  line  =  .5) 

hist(PR0Bt[,4],  col=rgb(0,0,l,.25),  axes  =  FALSE,main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 

mtext(side  =  4,  "Tri-Weibull",  line  =  .5) 

mtext(  expression("AIC"[corrected]*"  (40, 100, Low)" ),  outer  =  TRUE) 
dev.off(2) 

# 

# 

tiff(filename  =  "50byl00.tiff", 

width  =  1800,  height  =  2700,  units  =  "px",  pointsize  =  12, 
compression  =  c("none"),bg  =  "white",  res  =  300) 

par(mfrow  =  c(3,l),oma=rep(4,  4),  mar=rep(l,  4) ) 

hist(PROBs[,5],  col=rgb(l,0,0,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 

mtext(side  =  4,  "Weibull",  line  =  .5) 

hist(PROBb[,5],  col=rgb(0,l,0,.25),axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 

xaxs="r",  yaxs="r",  xlim=c(0,l),  breaks=g) 

axis(side  =  1,  tck  =  -.015,  labels  =  NA) 

axis(side  =  2,  tck  =  -.015,  labels  =  NA) 

axis(side  =  1,  lwd  =  0,  line  =  -.6) 

axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 

mtext(side  =  1, line  =  2.0) 

mtext(side  =  2,  "frequency",  line  =  2.5) 

mtext(side  =  4,  "Bi-Weibull",  line  =  .5) 

hist(PROBt[,5],  col=rgb(0,0,l,.25),  axes  =  FALSE, main="",  xlab  =  NA,  ylab  =  NA, 
xaxs="r",  yaxs="r",  xlim=c(0,l),  breaks=g) 
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axis(side  =  1,  tck  =  -.015,  labels  =  NA) 
axis(side  =  2,  tck  =  -.015,  labels  =  NA) 
axis(side  =  1,  lwd  =  0,  line  =  -.6) 
axis(side  =  2,  lwd  =  0,  line  =  -.6,  las  =  1) 
mtext(side  =  1, line  =  2.0) 
mtext(side  =  2,  "frequency",  line  =  2.5) 
mtext(side  =  4,  "Tri-Weibull",  line  =  .5) 

mtext(  expression("AIC"[corrected]*"  (50, 100, Low)" ),  outer  =  TRUE) 
dev.off(2) 

# 

# 
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